{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# _*Pricing Long Butterfly Options*_ \n",
    "\n",
    "The latest version of this notebook is available on https://github.com/Qiskit/qiskit-iqx-tutorials.\n",
    "\n",
    "***\n",
    "### Contributors\n",
    "Stefan Woerner<sup>[1]</sup>, Daniel Egger<sup>[1]</sup>\n",
    "### Affliation\n",
    "- <sup>[1]</sup>IBMQ"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Introduction\n",
    "<br>\n",
    "Suppose a <a href=\"http://www.theoptionsguide.com/butterfly-spread.aspx\">long butterfly option</a> with strike prices $K_1 < K_2 < K_3$, with $K_2 - K_1 = K_3 - K_2$ and an underlying asset whose spot price at maturity $S_T$ follows a given random distribution.\n",
    "The corresponding payoff function is defined as:\n",
    "<br>\n",
    "<br>\n",
    "$$ F(S_T) = \n",
    "\\begin{cases}\n",
    "0         ,& S_T < K_1 \\\\\n",
    "S_T - K_1 ,& K_1 \\leq S_T < K_2 \\\\\n",
    "2K_2 - K_1 - S_T ,& K_2 \\leq S_T < K_3 \\\\\n",
    "0         ,& S_T \\geq K_3. \n",
    "\\end{cases}$$\n",
    "<br>\n",
    "In the following, a quantum algorithm based on amplitude estimation is used to estimate the expected payoff, i.e., the fair price before discounting, for the option:\n",
    "<br>\n",
    "<br>\n",
    "$$\\mathbb{E}\\left[ F(S_T) \\right].$$\n",
    "<br>\n",
    "The approximation of the objective function is explained in detail in the following paper:<br>\n",
    "<a href=\"https://arxiv.org/abs/1806.06893\">Quantum Risk Analysis. Woerner, Egger. 2018.</a>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "import numpy as np\n",
    "\n",
    "from qiskit import BasicAer\n",
    "from qiskit.aqua.algorithms import AmplitudeEstimation\n",
    "from qiskit.aqua.components.uncertainty_models import LogNormalDistribution\n",
    "from qiskit.aqua.components.uncertainty_problems import UnivariateProblem\n",
    "from qiskit.aqua.components.uncertainty_problems import UnivariatePiecewiseLinearObjective as PwlObjective"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Uncertainty Model\n",
    "\n",
    "We construct a circuit factory to load a log-normal random distribution into a quantum state.\n",
    "The distribution is truncated to a given interval $[low, high]$ and discretized using $2^n$ grid points, where $n$ denotes the number of qubits used.\n",
    "The unitary operator corresponding to the circuit factory implements the following: \n",
    "$$\\big|0\\rangle_{n} \\mapsto \\big|\\psi\\rangle_{n} = \\sum_{i=0}^{2^n-1} \\sqrt{p_i}\\big|i\\rangle_{n},$$\n",
    "where $p_i$ denote the probabilities corresponding to the truncated and discretized distribution and where $i$ is mapped to the right interval using the affine map:\n",
    "$$ \\{0, \\ldots, 2^n-1\\} \\ni i \\mapsto \\frac{high - low}{2^n - 1} * i + low \\in [low, high].$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# number of qubits to represent the uncertainty\n",
    "num_uncertainty_qubits = 3\n",
    "\n",
    "# parameters for considered random distribution\n",
    "S = 2.0 # initial spot price\n",
    "vol = 0.4 # volatility of 40%\n",
    "r = 0.05 # annual interest rate of 4%\n",
    "T = 40 / 365 # 40 days to maturity\n",
    "\n",
    "# resulting parameters for log-normal distribution\n",
    "mu = ((r - 0.5 * vol**2) * T + np.log(S))\n",
    "sigma = vol * np.sqrt(T)\n",
    "mean = np.exp(mu + sigma**2/2)\n",
    "variance = (np.exp(sigma**2) - 1) * np.exp(2*mu + sigma**2)\n",
    "stddev = np.sqrt(variance)\n",
    "\n",
    "# lowest and highest value considered for the spot price; in between, an equidistant discretization is considered.\n",
    "low  = np.maximum(0, mean - 3*stddev)\n",
    "high = mean + 3*stddev\n",
    "\n",
    "# construct circuit factory for uncertainty model\n",
    "uncertainty_model = LogNormalDistribution(num_uncertainty_qubits, mu=mu, sigma=sigma, low=low, high=high)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZwAAAEyCAYAAADOV2anAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3debgcVZnH8e9PEAgEMLKERSSACiJRMEGMoCSArDpBEKKiYxCNjAKOIgwiQgAdAWVxYBQjakRGoyKiIIhhSRDZ40JYomwhBhQEAjGEQELe+eOcC51K973dt/tWNcnv8zz93FtVp6re7tu33646myICMzOzgfaKqgMwM7OVgxOOmZmVwgnHzMxK4YRjZmalcMIxM7NSOOGYmVkpKk84kraVdI2khZIekXSKpFX62OdNkn6Tyz8naY6kCyRtXCg3WVLUeWwzsM/KzMyKVq3y5JKGAFcDdwNjga2AM0mJ8IRedl0XeBC4EHgE2AI4CRghaceIWFJTdhZwaGH/2c3Et/7668ewYcOaKTognnnmGdZaa63Kzt9It8YFjq0/ujUucGz9UXVcM2bMeDwiNqi7MSIqewBfAOYB69SsOxZYWLuuyWO9GwjgrTXrJgO39ze+ESNGRJWuu+66Ss/fSLfGFeHY+qNb44pwbP1RdVy9feZWfUttH+CqiJhfs24KMAjYtcVjPZF/rtaJwMzMrLOqTjjbkG55vSgi5pCucPqsZ5H0CkmrSdoaOA24Dbi1UGxbSfNzXc8NklpNZGZm1gGKCsdSk7QYOCYizimsnwtcGBHH97H/b4C98uIMYN+IeKxm+2eA50l1RBsARwMjgF0iopiYevaZAEwAGDp06IgpU6b056l1xIIFCxg8eHBl52+kW+MCx9Yf3RoXOLb+qDquMWPGzIiIkXU3NrrXVsYDWAz8Z531c4H/bmL/1wM7AR8mXSnNANbopfyapMYGlzYTn+tw6uvWuCIcW390a1wRjq0/qo6LLq7DmUdqcVY0JG/rVUTcGxG3RMRFpCudHYAP9VJ+IXAF8Nb+hWtmZv1VdcKZRaGuRtJmpCuRWXX3aCAiHgKeBLbsq2h+mJlZiapOOFcCe0lau2bdOOBZYHorB8oNB9Yj3TJrVGYQsB/p1puZmZWo0o6fwPnAUcAlkk4nXZ1MBM6KmqbSku4DpkfEYXn568AS4BbgKeCNpP4795OaVSNpXeBy4CLgPmB94LPAJsBBJTw3MzOrUWnCiYh5knYHzgMuIyWPs0lJp9aqQO1wN7cDR5Jak60BzAF+Dnw1Ip7JZZ4D/kkasWBDYBFwE7BrRNw+EM/HzMwaq/oKh4i4G9itjzLDCstTyFcyveyzCDig3fjMAIYd9+u2j3H08CWMb/M4s0/br+04zKpSdR2OmZmtJJxwzMysFE44ZmZWCiccMzMrhROOmZmVwgnHzMxK4YRjZmalcMIxM7NSOOGYmVkpnHDMzKwUTjhmZlYKJxwzMyuFE46ZmZXCCcfMzErhhGNmZqVwwjEzs1I44ZiZWSmccMzMrBSVJxxJ20q6RtJCSY9IOkXSKn3s8yZJv8nln5M0R9IFkjauU3aspJmSFkm6W9K4gXs2ZmbWyKpVnlzSEOBq4G5gLLAVcCYpEZ7Qy67rAg8CFwKPAFsAJwEjJO0YEUvy8XcBfg58EzgK2Bf4saR5EfHbAXlSZmZWV6UJBzgcGAQcEBHzgamS1gEmSjojr1tORNwI3FizapqkucBvgTcDf8jrvwRcHxFH5eXrJL0JODGXNTOzklR9S20f4KpCYplCSkK7tnisJ/LP1QAkrQ6MAX5aKDcFGCVp3dbDNTOz/qo64WwDzKpdERFzgIV5W68kvULSapK2Bk4DbgNuzZu3Al5ZPD5wD+l5v6G90M3MrBWKiOpOLi0GjomIcwrr5wIXRsTxfez/G2CvvDgD2DciHsvbdgZuAHaIiD/V7PM64F5gr3r1OJImABMAhg4dOmLKlCn9fXptW7BgAYMHD67s/I10a1wwcLHNfPjpto8xdBA8+mx7xxi+aecvzFfGv2cndGtsVcc1ZsyYGRExst62qutw2nUk8Grg9aRGBldK2jkiFvX3gBExCZgEMHLkyBg9enQn4uyXadOmUeX5G+nWuGDgYht/3K/bPsbRw5dw5sz2/uVmHzK67TiKVsa/Zyd0a2zdGhdUn3DmkVqcFQ3J23oVEffmX2+R9DtSy7UPAd+r2b94/CE15zYzs5JUXYczi0JdjaTNgDVZvu6lVxHxEPAksGVedT+wuHj8vLwU+Gs/4jUzs36qOuFcCewlae2adeOAZ4HprRwoNxxYj3SVQ0Q8B1wHHFQoOg64KSLavylvZmZNq/qW2vmkDpmXSDqddHUyETirtqm0pPuA6RFxWF7+OrAEuAV4CngjcCzpqqa2lv9UUh+dc4BLSR0/9wX2HtinZWZmRZUmnIiYJ2l34DzgMlLyOJuUdGqtCtQOd3M7qcHABGANYA5pRIGvRsQzNce/QdL7gS8D/0Gu4/EoA7aiGNahxgztNoqYfdp+bcdhK76qr3CIiLuB3fooM6ywPIVlr2R62/dS0tWNmZlVqOo6HDMzW0k44ZiZWSmccMzMrBROOGZmVgonHDMzK4UTjpmZlcIJx8zMSuGEY2ZmpXDCMTOzUjjhmJlZKZxwzMysFE44ZmZWCiccMzMrRcujRUsaDrwN2Ig0NcCTpNkzb4wIT9tsZmZ1NZVwJG1Jmk/mEGAoaYrmp4DngFeRpoReKmk6cAHwk4hYOiARm5nZy1Kft9QkXQDcBWwPnALsAKwRERtExGsiYjCwIfBeYCZwBnCPpF0GLmwzM3u5aeYK51lgm4h4qFGBiHgcuBK4UtLngIOATTsTopmZrQj6vMKJiCN7SzZ1yi+NiJ9ExE+aKS9pW0nXSFoo6RFJp0hapY99dpT0fUn35f3+IukkSWsUyk2UFHUeezf7fMzMrDPammJa0nbAroCA6RExs8X9hwBXA3cDY4GtgDNJifCEXnYdl8ueDtwLvBk4Nf88sFD2aaCYYO5pJU4zM2tfvxOOpP8AvgJcA6wFfE3S0RHxzRYOczgwCDggIuYDUyWtA0yUdEZeV89p+TZej2mSFgHflrR54YpsSUTc3EJMZmY2AJppNLBmg03/BYyKiIMiYl/gCOCLLZ5/H+CqQmKZQkpCuzbaqZBsevwx/9ykxRjMzKwEzXT8/KukQ+qsF6l5dI/+NIPeBphVuyIi5gAL87ZWjMox3F9Y/ypJj0taLOmPkg7oR5xmZtYmRUTvBaR3AecAzwNHRcStef2nSc2kryH1w9kdODYizm365NJi4JiIOKewfi5wYUQc3+RxNgLuAK6IiPE16z9MarL9R2Bt4JPAvsCBEXFJg2NNACYADB06dMSUKVOafTodt2DBAgYPHlzZ+Rvp1rhg4GKb+fDTbR9j6CB49Nn2jjF803WXWe7WuDplZXyvtavquMaMGTMjIkbW29ZnwgGQJOAwUsX8VOC/IuLvkt7CS7e+ro+IP7USWCcSjqTVSA0PXgOM6G20g/w8bgQGRcT2fR175MiRcfvtt/dVbMBMmzaN0aNHV3b+Rro1Lhi42IYd9+u2j3H08CWcObOtdjrMPm2/ZZa7Na5OWRnfa+2qOi5JDRNOU2OpRXIBsDXwKHCnpC8CsyLif/KjpWSTzQPqfTUakrf1KieQC4E3Afv2NbROpOx6CfDmvppem5lZZ7U0eGdEzI+IY4CdSOOpzZL0/jbOP4tCXY2kzUi36GbV3WNZ55CaU4+NiGbKA0R+mJlZiZpqpSbpy5JuyZXuk4BFETGWVNdxkqTp+fZaq64E9pK0ds26caTRDab3EdcXSC3jPhwRNzRzsnxFdCDw54h4oR/xmplZPzVz4/a7wLakPjcLSUlmqqRtI2JqTjSfyusujYgJLZz/fOAo4BJJpwNbAhOBs2qbSku6j9Sx9LC8/CHgv4HJwMOS3l5zzPsj4p+53HTg56SrpbWAT5CuzvZvIUYzM+uAZhLOPsBBETEVQNLvgSdIPf3vy6NCnyfpR8BJrZw8IuZJ2h04D7iMNAL12aSkU4yzts5lz/xzfH7UOpSUiADuA/4T2JjUZPoPwH4RcWUrcZqZWfuaSTizgI9ImgEsIjUtfgaYW1soIp4EPtNqABFxN7BbH2WGFZbHs3yiqbffYa3GY2ZmA6OZhPNR0hXD46TK9tmkK55FAxeWmZmtaPpMOBHxF2CUpLWA1Tyrp5mZ9UfTvb0i4hnSrTQzM7OWNdMs+iOtdpKU9DpJ7+x/WGZmtqJppuPn54D7JZ3aW18bSetJOkTSZcCfSC3DzMzMgObqcHaQNA44EviipAWkCcweB54DXgVsAbyWNBzNRcDhEfHwgEVtZmYvO03V4eTpon8iaStgD+CtwEakzpSPAtcDvwemRcTiAYrVzMxexloaIjYi7mf5+WbMzMz61NLgnWZmZv3lhGNmZqVwwjEzs1I44ZiZWSlaSjiS3ivJScrMzFrWavK4FJgr6XRJbxyIgMzMbMXUasLZCvgOcDBwp6SbJH1C0jqdD83MzFYkLSWciJgdESdFxBbAu0kTnJ0N/F3SDyWNGYggzczs5a/f9TERcW1EfAR4AzADOAS4WtIDkj4rqaVOpWZmtmLrd8KRtKukycBfgO2A/yVN/XwxcDJwYScCNDOzFUOrrdQ2l3SipPuBa4HNgAnAxhFxZERcExHHkmYJHdvkMbeVdI2khZIekXRKX9MhSNpR0vcl3Zf3+4ukkyStUafszpJukbRI0oOSjmrlOZuZWWe0etvrAeAR0pTT34uIBxuUuwu4ta+DSRoCXA3cTUpQWwFnkhLhCb3sOi6XPR24F3gzcGr+eWDN8V8HXAVcDnwBeBtwlqSFEXFBX/GZmVnntJpw3gNcFRFLeysUEX8FmmlAcDgwCDggIuYDU3OLt4mSzsjr6jktIh6vWZ4maRHwbUmbR8RDef0xpAT54YhYAlwr6bXASZK+GxHRRIxmZtYBrdbh7EialmA5kjaWdGKLx9uHlMBqE8sUUhLatdFOhWTT44/55yaF41+Sk03t8V9DqncyM7OStJpwTiJ9WNezSd7eim2AWbUrImIOsDBva8UoYCl5+gRJa5HqmGYVyt1Tc24zMyuJWrmrJGkpsFNE3FZn21jguxGxfgvHWwwcExHnFNbPBS6MiOObPM5GwB3AFRExPq/bFJgLvC8iLq0puyqwGPhkREyqc6wJpIYQDB06dMSUKVOafTodt2DBAgYPHlzZ+Rvp1rhg4GKb+fDTbR9j6CB49Nn2jjF803WXWe7WuDplZXyvtavquMaMGTMjIkbW29ZnHY6kj5JanQEE8C1JxbqVNYDhwG/bCbQ/JK0G/BRYAHy23ePlJDQJYOTIkTF69Oh2D9lv06ZNo8rzN9KtccHAxTb+uF+3fYyjhy/hzJntdU+bfcjoZZa7Na5OWRnfa+3q1riguUYDC4En8u8CngaeLJR5HrgS+GaL558H1PtqNCRv65Ukkfr7vAnYOSJq93kq/ywef0jNuc3MrCR9JpyI+BnwMwBJ3wdO6aU5dKtmUahLkbQZsCbL173Ucw6pOfW7I6JYF/SMpL8Vj1+z3MzxzcysQ1odS+3QDiYbSFdFe0lau2bdOOBZYHpvO0r6AnAEqcnzDb0c/32FjqTjgL8Bd/Y7ajMza1nVc9ucDzwHXCJpj1xhPxE4q7apdB5R4Ls1yx8C/pt0O+1hSW+veWxQc/yvkVrV/VDSGEnHAp8kXaW5D46ZWYmaaTRwKzA+Iu6WdBup4UBDEfG2Zk8eEfMk7Q6cB1xGqnc5m5R0inHWXqXsmX+Oz49ah5JGQiAi7pO0N3AW6WrnH8DRHmXAzKx8zTQauIt0i6vn945eGUTE3cBufZQZVlgez/KJptG+N5CGtDEzswo102jg0Jrfxw9oNGZmtsKqug7HzMxWEs3U4fRZb1OrlTocMzNbeTRbh+MWXWZm1pZm6nDGlxCHmZmt4FyHY2Zmpai0H46Zma08Ku+HY2ZmKwf3wzEzs1K0PAlGnn9mPKn3/sbA34FbgB9ExPMdjc7MzFYYLTUakPRG4F7gf4HtgBfyz/8F7pO0bccjNDOzFUKrVziTSBOwvTMi5vSslPRa4HLS6M/v6lx4Zma2omg14YwEPlibbAAiYo6kk4AfdSwyW+kM69B0ye1Ouzz7tP3ajsPMltdqP5zZwBoNtq0BzGmwzczMVnKtJpzjgC9L2ql2paS3A6cC/9WpwMzMbMXSn8E71wFulPQY8BiwYX48ARwPXDoAcZqZ2ctcfwbvvGuAYjEzsxVY5YN35qbU5wKjSFNMXwCcHBEv9LLPasBXgLeTGjKsERGqU24y8NE6h3hjRMxqP3ozM2tWyx0/O0nSEOBq4G5gLLAVcCapbumEXnZdE/g4cCtwI71PUT0LOLSwbnb/IjYzs/6qNOEAhwODgAMiYj4wVdI6wERJZ+R1y4mIpyS9OiJC0hH0nnCeiYibOx+6mZm1ouXpCSSNk3S1pDmSHis+WjzcPsBVhcQyhZSEdu1tx4jwIKJmZi8jrQ5t8yHgB8B9wGuAX5FGGHgFMB84r8Xzb0O65fWi3Kl0Yd7WCdtKmi/pOUk3SOo1kZmZ2cBQKxcKkv4IXAycBiwGRkbEHyStDUwFLo6Ir7dwvMXAMRFxTmH9XODCiDi+iWMcAZzboNHAZ4DnSXVEGwBHAyOAXSLi1gbHmwBMABg6dOiIKVOmNPt0Om7BggUMHjy4svM3MlBxzXz46baPMXQQPPps3+V6M3zTdZdb162xdWtcndKt/wPQvbFVHdeYMWNmRMTIettarcN5PfD7iHhB0gukPjlExL8knQ6cDTSdcAZaRHyjdlnSFaRm3ccD+zfYZxJpzDhGjhwZo0ePHuAoG5s2bRpVnr+RgYqr3SFpIA1tc+bM9qomZx8yerl13Rpbt8bVKd36PwDdG1u3xgWt1+HMB1bPvz8MvLFmm4D1WjzePKDeV6MheVtHRcRC4ArgrZ0+tpmZ9a7VrzW3AW8GriLV35woaQnpttWJQKutwWZRqKuRtBmp2fNA9ZMJPGupmVnpWk04XwU2z7+fmH//FulK6Tbgky0e70rgGElrR8S/8rpxpCmtp7d4rD5JGgTsB8zo9LHNzKx3LSWc3J/l5vz7U8BYSasDqzfqM9OH84GjgEtyHdCWwETgrNrjSboPmB4Rh9Ws2wdYC9g+L78/b7otIh6StC6pBd1FpFZ16wOfBTYBDupHrGZm1oaOTTEtqeUppiNinqTdSc2pLyMNbXM2KekU41ylsO5bvHS1BfCz/PNQYDLwHPBP0ogFGwKLgJuAXSPi9lbiNDOz9rWUcPIU078hXSXMII0WvR3w78CXJO0dEXe3csxcvreRAoiIYc2sK2xfBBzQSixmZjZwPMW0mZmVotVm0SOBE+tNMQ2cBOzYqcDMzGzF4immzcysFK3eUjsOOFPSgxFxS8/KmimmP9/J4Mzs5WtYh0ZBaHc0hdmn7dd2HNYZnmLazMxK4SmmzcysFJVPMW1mZiuHfg0RK2kTYBTwatKttJsj4pFOBmZmZiuWVjt+rgKcC3yCZXv+vyBpEnBkRCztYHxmZraCaLVZ9MnAx0iNA4aRpoIelpc/xvJD0piZmQGt31L7d+CEwqyec4CvSQrSQJwndio4MzNbcbR6hbMhcEeDbXfk7WZmZstpNeH8FfhAg20fAP7SXjhmZraiavWW2peBKXmwzouBR0lXNQcBY2icjMzMbCXX6gRsP5X0FKnxwDeAVwKLSVMV7B0RUzsfopmZrQiaTjiSXkmadO3OiBgl6RWkWTQfd1NoMzPrSyt1OC8A1wLbAETE0oh4zMnGzMya0XTCyYnlXmCjgQvHzMxWVK22UvsicKKk4Z0KQNK2kq6RtFDSI5JOySMa9LbPapK+Jul3kp7NfYAalR0raaakRZLuljSuU7GbmVnzWm2ldgKwHvAnSQ+TWqkt82EfEW9r9mCShgBXA3cDY4GtgDNJifCEXnZdE/g4cCtwI7Bbg+PvAvwc+CapU+q+wI8lzYuI3zYbp5mZta/VhHMXcGcHz384aXicAyJiPjBV0jrAREln5HXLiYinJL06IkLSETRIOMCXgOsj4qi8fJ2kN5FGQ3DCMTMrUavNosd3+Pz7AFcVEssU4HRgV+CyXmJpeBsNQNLqpL5BRxU2TQG+L2ndiHi6X1GbmVnLmqrDkTRI0oGSjpb0IUlDO3T+bYBZtSsiYg6wMG9rx1akfkKzCuvvIT3vN7R5fDMza4H6uFBA0pakepZhNavnAwe3Ww8iaTFwTEScU1g/F7gwIo5v4hhHAOdGhArrdwZuAHaIiD/VrH8dqbXdXvXilzQBmAAwdOjQEVOmTGn9iXXIggULGDx4cGXnb2Sg4pr5cPsXnEMHwaPPtneM4Zuuu9y6bo2tW+OC7o6tE1a2/89mjRkzZkZEjKy3rZlbamcAS4F3kkYU2IJUCf/t/PsKJSImAZMARo4cGaNHj64slmnTplHl+RsZqLjGH/frto9x9PAlnDmzX/MKvmj2IaOXW9etsXVrXNDdsXXCyvb/2QnN3FIbRZqS4PcRsSgi7gE+CbxW0sZtnn8eUO/rx5C8rd1jU+f4QwrbzcysBM0knI2BBwrr7gdE+51AZ1Goq5G0GanZc7HupVX3k8Z5K9YFbUO6Yvtrm8c3M7MWNNvxs/eKnv67EthL0to168YBzwLT2zlwRDwHXEcaybrWOOAmt1AzMytXszdHr5K0pM76a4rrI6KVSdjOJzVbvkTS6cCWpGmqz6ptKi3pPmB6RBxWs24fYC1g+7z8/rzptoh4KP9+KjBN0jnApaSOn/sCe7cQo5mZdUAzCefkgTp5RMyTtDtwHqnPzVPA2aSkU2tVoDjczbeAzWuWf5Z/HgpMzse/ISeiLwP/ATwIfMijDJiZla/PhBMRA5Zw8vHvpvFIAT1lhjWzrsG+l5KubszMrEKtDt5pZmbWL044ZmZWCiccMzMrhROOmZmVwgnHzMxK4YRjZmalcMIxM7NSOOGYmVkpnHDMzKwUTjhmZlYKJxwzMyuFE46ZmZXCCcfMzErhhGNmZqVwwjEzs1I44ZiZWSmccMzMrBSVJxxJ20q6RtJCSY9IOkVScTrpevutK+n7kuZJelrS/0lar1BmsqSo89hm4J6RmZnV0+cU0wNJ0hDgauBuYCywFXAmKRGe0MfuPwXeAHwcWAqcTppK+p2FcrOAQwvrZrcTt5mZta7ShAMcDgwCDoiI+cBUSesAEyWdkdctR9IoYE9g14i4Pq97GLhF0h4RcXVN8Wci4uaBfRpmZtaXqm+p7QNcVUgsU0hJaNc+9nu0J9kARMStwIN5m5mZdZmqE842pFteL4qIOcDCvK3p/bJ76uy3raT5kp6TdIOk3hKZmZkNEEVEdSeXFgPHRMQ5hfVzgQsj4vgG+00l3Srbv7D+ImDLiHhHXv4M8DypjmgD4GhgBLBLviKqd+wJwASAoUOHjpgyZUobz7A9CxYsYPDgwZWdv5GBimvmw0+3fYyhg+DRZ9s7xvBN111uXbfG1q1xQXfH1gkr2/9ns8aMGTMjIkbW21Z1Hc6Aiohv1C5LugK4Czge2L/BPpOASQAjR46M0aNHD3CUjU2bNo0qz9/IQMU1/rhft32Mo4cv4cyZ7b2tZx8yerl13Rpbt8YF3R1bJ6xs/5+dUPUttXlAva8fQ/K2ju4XEQuBK4C3thCjmZl1QNUJZxaFOhdJmwFrUr+OpuF+WaO6nVqRH2ZmVqKqE86VwF6S1q5ZNw54Fpjex34bSdqlZ4WkkcCWeVtdkgYB+wEz2gnazMxaV3XCOR94DrhE0h65wn4icFZtU2lJ90n6bs9yRNwE/Ba4UNIBkvYH/g+4oacPTh6J4HeSPilpd0njgOuATYD/LusJmplZUmmjgYiYJ2l34DzgMuAp4GxS0qm1KlAc7mZcLvs9UuK8HDiqZvtzwD9JIxZsCCwCbiJ1Fr29o0/EzMz6VHkrtYi4G9itjzLD6qx7ijRkTXHYmp7ti4ADOhCima1ghnWoBV07LfFmn7Zf2zG83FR9S83MzFYSTjhmZlYKJxwzMyuFE46ZmZXCCcfMzErhhGNmZqVwwjEzs1I44ZiZWSmccMzMrBSVjzRg5eqGHtawcvayNlvZ+QrHzMxK4YRjZmalcMIxM7NSOOGYmVkpnHDMzKwUTjhmZlYKJxwzMyuFE46ZmZWi8o6fkrYFzgVGAU8BFwAnR8QLfey3LnAOsD8pcV4OHBURTxTKjQW+DLweeCAf+yedfh5mZu1a0TtmV3qFI2kIcDUQwFjgFOBo4OQmdv8pMBr4ODAe2BG4tHD8XYCfA9cB+wC/Bn4sac+OPAEzM2ta1Vc4hwODgAMiYj4wVdI6wERJZ+R1y5E0CtgT2DUirs/rHgZukbRHRFydi34JuD4ijsrL10l6E3Ai8NuBe1pmZlZUdR3OPsBVhcQyhZSEdu1jv0d7kg1ARNwKPJi3IWl1YAzpSqjWFGBUviVnZmYlqTrhbAPMql0REXOAhXlb0/tl99TstxXwyjrl7iE97zf0I14zM+snRUR1J5cWA8dExDmF9XOBCyPi+Ab7TQWeiYj9C+svAraMiHdI2hm4AdghIv5UU+Z1wL3AXhGx3G01SROACXlxa+Av/X6C7VsfeLzC8zfSrXGBY+uPbo0LHFt/VB3X5hGxQb0NVdfhdJ2ImARMqjoOAEm3R8TIquMo6ta4wLH1R7fGBY6tP7o1Lqj+lto8oF5dypC8rZ39en4Wyw0pbDczsxJUnXBmUairkbQZsCb162ga7pfV1u3cDyyuU24bYCnw137Ea2Zm/VR1wrkS2EvS2jXrxgHPAtP72G+j3M8GAEkjgS3zNiLiOVL/m4MK+44DboqIp9sPf8B1xa29Oro1LnBs/dGtcYFj649ujavyRgNDgLuBO4HTSQnjLOCciDihptx9wPSIOKxm3VWk0QM+T7piOR14LCLeWVNmF2AacB6pU+i+ufze9RoMmJnZwKn0Cici5gG7A6sAl5FGGDgbOKlQdNVcptY40lXQ94ALgRnA+wrHvwF4P7AHcBXwb8CHnGzMzMpX6RWOmZmtPKquwzEzs5WEE1I1XdsAABkPSURBVI6ZmZXCCcfMzErhkQa6gCSRGjzsB7wReDWp5d0/gJuByRFRSb+h3C9qX0DAzyLiCUmvIbX22wqYDUyKiJklxvRfwBVlnrNZkgYBq0bEv2rWbQAcAWxL+rv+Cfjmy6RpfiXy/8R7gbeSpi+5nfQ374pK5zyq/ePAbrlxUlUx7AasBvw6Ip7J77VPk1r8PkD633ykivjqcaOBiuU3yBXACFKCeR7YlPRPdiXpjbM1cGpEnFpybG8DpgJrAUuAJ4G9crwvAHcB2wEbAXtExO9Kimsp6fWZBfwImBIR95dx7r5IugK4NyI+k5dHkf6OS0ktKUX6Wz9P+rC6q6S4dgAGRcSNNev2Br7AS4nwz8DE2jIlxXYjcFhE3JOXh5CmDxkBLMjFBpO+fO1Vm8wHOK5P9bJ5EPA14BuksRmJiG+WERe8OCbkNcBmedWDpClbpgKvInV835rUp3FERMwtK7ZeRYQfFT6AH5PesMNr1m0C/Ab4eV7elfSP97GSY5tK6jz7KtLI2+cBc4FfAq/MZVYnfaBeV2JcS4GvkmZ5fY6U/G4DPgtsWvHf83FgbM3yzaQPhrVr1q1LatJ/VYlx3Qx8sWb5Y/l1vAb4InBC/lsvqY2/xL/n22qWv0v6crN3zbq9ScNRnV1yXC/kn/UetdteKPk1+ynpC8LrSHdEfpg/R27sea+RBvH8M/DtMmPrNe6qA1jZH6RptQ+ss35YfkNvnJePB/5ccmxPAPvULG+Y/7n2LJTbD3i8xLhe/IDK/2wT8gfnkvyYlte9uoK/50LgXTXLzxdfr5rX7JkS45pfGwdwH3BunXLnV/A+KyacfwL/Wafc54GHSozrF8DfgUPJd4Nqtr0qx/2usuIpnP8R4OCa5c1zPAcUyh0K/LWKGOs93Gigeq8gJZaiF0i3X3oGH72FaubwiTq/F+/DVnZfNiKejIhJEbE78BrSFOWrkT44/y6p/UniW3MnaeK/Ho+SkmLReqTkVJalheXNgYvrlLuYdCumSq8i1dkUzSDdvi1FRLwP+ChwDHBbnvLkxc1lxdHAENIt+B4P558PFco9QPq/6ApOONWbCnxZ0pY9K/I97P8hvaF6GgsMBsquZL4d+LykwZJeQbrKehj4D0mr5FhXBT5F+qCtVET8IyK+ERHvALYgjVixSclhnAYcJ+lj+bX5CvA1Se+WtJqk1XPdyVdZfjbagfQ74JCa5buAekPY78hLH15lOlDSp3K9yTyg3nwq65Ou1EoTaVSSN5NmCv61pCm50UzVHiN9aejxAvBt0hecWhsCpdR5NaXqS6yV/UH69nEXaWTr+0hjyz1LutVWezvrDOAnJcc2kvTPvzjH9ATwFlISfIA0HNGDpHqUMSXGtcwtmG57AB8nfTA+Ddyaf3+BdLvvhfz4BbBmiTENz3H8EHgbaSr2x0gJ8d2kCufTgEXUuZ1Vwt+z+PhenXLfBn5X4d91I9IwWguAM/PfsapbapfWe43qlDsXmFrVa1Z8uJVaF8hXCweTPszXICWeH0XEk5UGBuRvc+8hNaH/eUT8XdJGwLGkWy8PARdExB9KjOkk4DvRRc09iyStRxrv722kDyqRkvc9wOURMaOCmLYHvgXsRLolpLyp5/d5wCkR8Y2yY2uGpE8A90fEtRXHMYo05uPWwH5Rcqu+HMNQ0heWB/so9zlSndw15UTWOyccs5WMpDeSkk4xEd4YEYurjM1WbE44XUTSm0gTxNXOSjorSuqr0SpJr4iIYmV0ZSStQeqMuhS4r+oPz1yHsyU1HXkjYk6VMb3c5A6gRIUfVLkzryJiYc267ckdn6u4Wn25cqOBLpArmB8C7gB+RppAaVL+/Q5JsyUdWlFsB0i6VNIVkt6b142TNBtYLOmhfKujzJg+LOljNcurSjqN1HfjDlIDhiclHVdmXDXxjJD0K1Jl7T3A74GbgAclPSzpFElrVhFbN5K0Z2ESRiTtL+kPpPrD5yXdLmm/kuNaV9IvSHVf8yV9R9Iqkn4A/IH0/3mrpBskrV9mbM2SdKCkeq1gK+GEUzFJR5IqQy8HRpNalbwyPzYkdfq8HDhf0qdLju1gUjPZ9Un/+D/JyeWHpH4vR5E6mp0vaa8SQzue1OG0x+k5lq8C7yK9ZmcCJ0k6vsS4kLQn6TXZhHSf/1RSr/kXgImkCQYPBG7MrRHLjO09kq6RdI+kX0p6V50yO1XwAXUlaUinnhjeB1xCasBwXH4sBn6ZX9+ynAq8E/gcqaPsO0gtC3cjdUQdSqrf3CKXtT74llrFJD0AnB8RZ/RR7ljg8IjYsrdynSTpNmBGRByelw8hTXh3XkQcXVPu+8BmEbFHSXEtJLXgm56XHwO+UqzslvR54MiI2LzOYQYqthnAnRHx0cL6I0l9hLYk9RO6Ebg5InobPqWTcb2bNHrFzcAfgVHA9sA5wOd7bllJ2olUl1Oc8HAgY1sKvD0ibs3LfwAejoj3FspdAawVEbuWFNds0vvqO3l5B1JfoEMj4gc15T4BHB8RW5QRVz7n95osujkwusy/Z298hVO9jUlNZ/tyKyV2esu2ZtnOgZeTrryKnSkvIY3HVZanSVddPdYlDeFR9GfSVWKZtgUuqrP+IuC1wNYRsYj0Qf++OuUGyknAhRGxc0QcEREjgE8AnwQuyfVf3WI70lV/0STSYJ5l2YCX+sFBHjONNE5Zrfuo329oIH2U1JR9eB+P0r5sNcMJp3p/Bj6RO1bWlStOP0GqnyhTsOzU3j0DKT5VKLeA1Du8LL8idUhdLS9fDXywTrkPkvo4lekxUvP2oreQXs+ezrsP8dIoEmXYjkIijIjvkW4/vh24VlK9ERHKUnur5Wnqd1Z8hnI/sx4kvT493klq/PGOQrmdgbIbg9wLXBsRO/b2IN2O7BqenqB6R5Nuddwt6RLSCMg9H+jrklqtvY/UQXTvkmN7iPSN/SqAiHgh90G4p1BuK9KYU2X5Aqnn/J2SLiB1QD1d0nakcdREGl5mB9IQ92WaBJwqaTApET5P6r3/RdIApz19h7ak3A+pRaRRv5cRETPykC1XkW7zTSwxplpXSVqSf1+X9LebXiizDeW+z84HviFpOCkJHkx6730p/33/TLri+izl1+HczPKJr57a/laVc8KpWET8PjexPJY09MhmhSJ/I1Wqfi3KH4L/EgrjMEXELXXKfQAobU6QiHhS0ttJH+Kf46XbZqPy43nSkEHvjIjbyoorx/aVXCdxHHBiz2rSqOD/WVN0MfDfJYZ2B2l0gV8VN0TEAznpXAFMLjGmHifXWfdYnXUHkka0LkVEnJfvPHyQ1DDg2Ig4X9Jc0tBTPePhnQ98vay4snNJLeX6Mp1lx/arlBsNdJncXLbn9tRTtW3/u5Wk15JiLXWcq5rzD2PZToz3d0EfnFeSrvzWAB6o6rWpieeTpNZ9OzQawULSWqQhd/aICN9u70W+zb1+RPyz6lheTpxwzMysFL6l1iWUpnLehPTt/PE629cH9o2IC0sPro58D/sPwCFl37ZSl0/jrC6clrvb6WUyXXK+sqmd+noGKd7Sv7lLGkm6zSjSNPSzJL2FdIuy5312XkRcVXZsjfgKp2KSVie1Hjogr1pKGpH2c7UflhX1j9i3l81rAT8h1VXcCRARV5QUV1dO45xj6cppuZulNM7aQRFxSonn7MrpktWlU1/nWPYiNZZ5ktR6bwNgLKne9R5SX6sRpAYrB0bEpWXF1qsyhqT2o9fhw08ktUr7BGk6gKNIc1rcC7y+ptxOlD+NbVdOsUuXTuOcz9uV03K3EP+BFbzPunK6ZLp06ut83t+ThtZZJS8fn+P4bqHcD0kdjCt/b0V4iunKH6Rm0EcU1m0EXE+aandUXldFwrmdl6bY3bzweHP+hzy4Z12JcXXlNM75nN06Lfdrm3wcXsH7rCunS66TcLpi6ut8zqdJV8g9y0NyvLsVyu1JatBTWmy9PVyHU73NKHTojIh/SNqd9O3k6jykTJn9D3rsSLryOp3U7+XzkeffkNTTafEfEVGc1nag9UzjfH1e7pZpnHt047Tcs5s8p5os10kvl+mSu2Lq6+xZlu1X1fP7oEK5NUl9sLqCE071HgFez0sfngBEatb7AUnfIF06l95YINJXpEmSfgp8GZgp6dz8e5VOA/5P0t9Ir0vPNM5PkG6j9XT8LHsaZ3hpWu4bSMmudlruayN1nq1iWu5/AdcCF/RRbhdSn7AydfN0yQfmynnooqmvSbfUTpR0bz7310mzBf+XpOsj4l/5S+GxpITYHaq+xFrZH6TBMKf1UeYLlFxP0iCON5N68j8CfIZqp9jtummcc1zdOi33VNJQKH2Vq6IOpyunS6aLp74m1XfNrnmv30+6JdrzvzCTlJznAduXGVtvD1/hVO+bwDhJr44GHfIi4qtK8+W8u9zQlovjDmC0pA8AZ1DhkBkRcUGeq6RnGucn6YJpnCPi9jwUSnFa7nfx0rTcV1LytNykK+gJTZT7J4Wr7RJ8knTrpy8PkpJTKaL5zq+3k1psliYi7stDOe1MapxyTUQ8K2k06cvY1qRb8j+Kklr1NcPNoq1f8m2htYAFEdE1EzyZWfdywjEzs1J4vKSXiTy97XerjqOebo2tW+OC7o6tW0m6WtI1VcdR1K1xQffF5jqcl48xdO8XhG6NrVvjgi6NTdLVpDsfu1cdSx2iC18zujcu6LLYfEvNzF6Uvw2/IiK6Zkh7W3F0Teaz3klaI08D0HW6NbZujQu6N7aI2L1bk42kV3bja9atcUH3xeaE8/KxH6lZaDfq1ti6NS7o0tiq+oCS9GlJ90t6VtKfJX2kTrG3UvJr1q1xdXtsjTjhmK0kuvUDKvfrOpc0COuXSJ0YJ0u6WNIaZcbycoir22PrjRsNVEzStU0WrTekxoDq1ti6NS7o3thqPqB+TBq6/h2kD6ixwIcjosrxtj4PfD0iXhxSJ48l+H/AdZLeExFPOK6XTWwNudFAxSQtAf5CGgepN5sCO0W58+F0ZWzdGhd0b2ySbicNbVPvA+pB4D2RJoqrYt6lfwHvjYhphfXDSKMyrEKaBmCDMmPr1ri6Pbbe+AqnencBsyJiXG+FJL2fkofPoHtj69a4oHtj25r0rfhFEXGNpLeTPqBukrR3ifHUepo0AOYyImK2pHcAvwZuAk51XC/q5tgach1O9W4G3t5EuaD8scu6NbZujQu6N7aGH1Ck22uPkz6gdiwxph4zgP3rbYiIecDupPHK/qfMoOjeuKC7Y2vICad6ZwBHNlHuCmCLAY6lqFtj69a4oHtj6+YPqIuALSXVm9OIiHgW+DfS1ApzHBfQ3bE15Docs5WApIOAz5LqauqOSi5pFeBbwLsjouxEbSsBJxwzMyuFb6mZmVkpnHDMzKwUTjhmZlYKJxwzMyuFE471StJ4STMk/UvSPEl/lHTWAJ3rYEnjmyg3UVLUPB6R9HNJWzV5nsm5533lmn3OuWzP8763wfZ78/aJAxVDi8dd5nXu9HkkvULSEfk9+ayk+ZLukvQ/kvrVx0nJnyR9tMH2ybk3f71t58mT6vXKCccakvQFUjv+q4ADgH8Hfklq3z8QDgbGN1n2aWBUfnwe2B64RtJaTex7agvnGWitPGeARcAWkkbWrpS0IzAsbx/oGJpVfJ07fZ6fAF8GLiG9Jz9K6t/0juh/89uDgVcDP+rHvl8HDpH0un6ee4XnoW2sN0cA346I42vWXSbp5KoCqrEkIm7Ov98saQ7wO2Bf4GfFwrmPySoR8XxE3F9inJ32DPAH4AOkjpo9PgBcC4yoIqgeZb3OkvYB3g/sGxFX1mz6RX+vbrKjgB9GxOKac61KSp4fATYBPijpfuDkiHhxeKI8rMwNwH8AR7cRwwrLVzjWm1cB/yiurP322HPbRNL+kmZJWiTpBknbFvfLt1RmSnpO0t8kfSX/MyNpMnAgsGvNrbKJLcQ6I/8cVieuu0jf/Heq3VaI7V2SrpO0QNLTkqZJ2qFm+zslTZe0UNITkr4jae3eApI0StKvJP1d0jP5Vs0hta9dP5/zFODgng/W/PPgvL5jMeTX4OLC8UbnMtvVvpZ9vc6NziNpX0lLJW1ROM8Wef3YBq/BrvnncqNz9/fqJl+ZvAO4uLDpM8CxpFEYrgA+BnwPWK/OYX5OusrxZ2sdvsKx3vwBODJfPVzey3DnmwNnkebleBY4GbhK0ut7hr2XtCfpFsiFwDHAm0nfGtcDDs+/v5aU5D6Vjzu3hViH5Z//KKw7Azglr687z4uk0cBU4DrSbZlngJ1JIzr/UdLOwNXApaRv1esBpwFD8nIjmwO/B84nfRDvDHxf0tKI+DH9f86XkEYE2IV0VfdO0qjAlwBfKymGWsPo+3VudJ6/A4+QXveJNeXHA4+RBqGs55n882uSzoyIh1qMuZ7d83H/XFi/K2mk7TPyF6nf5zHo6rkRGAoMr3Mciwg//Kj7ICWFB0gDTS4ljYR8CrBOTZnJefs7atZtDiwBDq9ZdzNwXeH4xwIvAK/JyxcD05qIayJpsMlV8+MNpGQxH9i4ENf2dfafDNxes3wT6faUGpzvd3Vi3y0ff7smX0vlWL9N+vDqWd/Uc6593vn3XwL/m3//JnBp/v1xYGInYgCmARcX1o2ufd4tvs6NzvNlUpJSTZyzSfO9NHotNgLuyOcO4E7geGBwG+/3ScBtddZ/G/hbPudkYFgvx1g1v/c/0d84VuSHL/usoYi4A3gjqUL2m6QPgi8Bt0saXFP0sYi4sWa/h0i3uN4GL97XfyvL1638hHRbd1Q/wlsPWJwffwG2BMZFxN9ryjwcEX/q7SC5kcFOwA8if2IUtq+Z4/uppFV7HsAN+dwN60wkDVFqMfVQTawTSAmyXVOA90tanXSVtdzttBJi6NHn69yH75G+pIzOy2Py8vcb7RAR/wB2APYiXe29CvgKcKOk1QAkfTTfQvxTvo07K/8+Q9Ir6xx2I1LCLvoK6crnQdL/wufzVW+9uJYAT+VjWYETjvUqIp6LiMsi4oiI2Bb4OPB64LCaYo/V2fUxYOP8+/rAK4FHC2V6luuOeNuHp0lD6Y8EXkP61nlloUzxfPUMISXSv/eyfRVSwl1c83iO9Jw26+XYk4FxpNtce+Z4vwd0YgrgXwGDSR+GawGXVRBDj2Ze54Yi4gHS1dShedWhwK0RcVcf+70QEb+NiE+Rbtd9n3Qra1Te/oOI2J70ZWcJsHNEbB8RI6KmUUCNNUh/1+J55uTjvo90xb8LcIMadw94js6+visM1+FYSyLiu5LOALapWb1hnaIbkm7BQfrWuLhOuaH5Z93Ri/uwJCL66kvTTOXxPNLtwo0bbH8qH2ciqcK46JF6OynNK/8e4NMRcX7N+o58yYuIZyRdThoB+mcR8UyxTAdiWASsVlg3pF44TR6vNxcA31Fqin8ALbbyioilkn5LSlbFD/vXA/Oi7ymXn6TBlUlOUL9Rmqp7Immqh7MlnZMTUq1X0b/39ArPVzjWkKTlEomkDYB1WfZb7YZKswz2lHkt6VvlrZC+iZJusR1UONzBpA/7m/Ly85T8zTB/UN8C/HtPq686228Gto6I2+s86iYcYHXS/9eL35hzq7ZiH6Z2nvO3SFc25zfY3m4Mc1n2iwWkq6T+6u25XpK3TyHFXPcWIYCkoQ02/RuwkPT3rPUWmqvA/wt15iiq974Abss/X10ouwGwJvDXJs630vEVjvVmpqRfAr8l3SLbnNTJciHwg5pyjwMXSTqBl1qpPUa6ndPjJFLLte+TPkyGk1oufScielpFzQLGStqf9GH3SC8f6J10HKkV2pWSJpHu148iVXhfTmrccI2kpaSK73+RbuHsB3wxIpb7cImIpyXdBpwoaT4psR5HuhW4Tk3Rfj/nSPPZT+tle7sx/AI4TNLZpNZiY4B2pqFu+FwjYpGk/wM+Dfw4Ip7q5Tg/lfQv4KekxgUbAocAY0mV9cV930JqYNCX35Neqw0i4p81638k6Y/A9aTblyNIV5YPA/cUjjGSdMV3I7a8qlst+NG9D9I//29Jt40Wkf65fwRsU1NmMqmF1wGkb3XPkf5xl2u9RapLmEn6JjuXVP+was329Ukfck+Sb2M1iGsiubVWL7FPpqaFVF/bSE1frycl06dIrd62r9m+E/AbUku4Z4C7SU3B1+0lhtcB1+Tyc0iJa5nYm33OLTzvZVqptRsD8AVSC61/kWaZ/DeWb6XW1Ovc13MF9sjr9+jjOX4s/y3m5vfSk6SEOLpB+cuADzTxfl8NeAL4SGH9+/L5/kFK2vNJiX6HOsf4BoUWjX689PAEbNaW3KFvu4gY2VdZs97kusGDgS0jYmkHjzsH2Csiilcj9cp+A3hdROzXYPtkUqKcXWfbKsBDwHERcVFbQa+gfEvNzColaWtgW9KQMCd3ONkMIXWKbbZO5WvAXyW9IercKu3DQaRbyg3rn1Z2bjRgZlX7NulW7RWk4WM6JiLmRcSgSA1Xmik/l3TLrlGrxUtJt1zrEXBYpL44VodvqZmZWSl8hWNmZqVwwjEzs1I44ZiZWSmccMzMrBROOGZmVgonHDMzK4UTjpmZleL/AW46+1xJ9j1lAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# plot probability distribution\n",
    "x = uncertainty_model.values\n",
    "y = uncertainty_model.probabilities\n",
    "plt.bar(x, y, width=0.2)\n",
    "plt.xticks(x, size=15, rotation=90)\n",
    "plt.yticks(size=15)\n",
    "plt.grid()\n",
    "plt.xlabel('Spot Price at Maturity $S_T$ (\\$)', size=15)\n",
    "plt.ylabel('Probability ($\\%$)', size=15)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Payoff Function\n",
    "\n",
    "The payoff function equals zero as long as the spot price at maturity $S_T$ is less than the strike price $K_1$, then increases linearly, and is bounded by $K_2$.\n",
    "The implementation uses two comparators, that flip an ancilla qubit each from $\\big|0\\rangle$ to $\\big|1\\rangle$ if $S_T \\geq K_1$ and $S_T \\leq K_2, and these ancillas are used to control the linear part of the payoff function.\n",
    "\n",
    "The linear part itself is then approximated as follows.\n",
    "We exploit the fact that $\\sin^2(y + \\pi/4) \\approx y + 1/2$ for small $|y|$.\n",
    "Thus, for a given approximation scaling factor $c_{approx} \\in [0, 1]$ and $x \\in [0, 1]$ we consider\n",
    "$$ \\sin^2( \\pi/2 * c_{approx} * ( x - 1/2 ) + \\pi/4) \\approx \\pi/2 * c_{approx} * ( x - 1/2 ) + 1/2 $$ for small $c_{approx}$.\n",
    "\n",
    "We can easily construct an operator that acts as \n",
    "$$\\big|x\\rangle \\big|0\\rangle \\mapsto \\big|x\\rangle \\left( \\cos(a*x+b) \\big|0\\rangle + \\sin(a*x+b) \\big|1\\rangle \\right),$$\n",
    "using controlled Y-rotations.\n",
    "\n",
    "Eventually, we are interested in the probability of measuring $\\big|1\\rangle$ in the last qubit, which corresponds to\n",
    "$\\sin^2(a*x+b)$.\n",
    "Together with the approximation above, this allows to approximate the values of interest.\n",
    "The smaller we choose $c_{approx}$, the better the approximation.\n",
    "However, since we are then estimating a property scaled by $c_{approx}$, the number of evaluation qubits $m$ needs to be adjusted accordingly.\n",
    "\n",
    "For more details on the approximation, we refer to:\n",
    "<a href=\"https://arxiv.org/abs/1806.06893\">Quantum Risk Analysis. Woerner, Egger. 2018.</a>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "# set the strike price (should be within the low and the high value of the uncertainty)\n",
    "strike_price_1 = 1.438\n",
    "strike_price_2 = 1.896\n",
    "strike_price_3 = 2*strike_price_2 - strike_price_1\n",
    "\n",
    "# set the approximation scaling for the payoff function\n",
    "c_approx = 0.25\n",
    "\n",
    "# setup piecewise linear objective fcuntion\n",
    "breakpoints = [uncertainty_model.low, strike_price_1, strike_price_2, strike_price_3]\n",
    "slopes = [0, 1, -1, 0]\n",
    "offsets = [0, 0, strike_price_2 - strike_price_1, 0]\n",
    "f_min = 0\n",
    "f_max = strike_price_2 - strike_price_1\n",
    "butterfly_objective = PwlObjective(\n",
    "    uncertainty_model.num_target_qubits, \n",
    "    uncertainty_model.low, \n",
    "    uncertainty_model.high,\n",
    "    breakpoints,\n",
    "    slopes,\n",
    "    offsets,\n",
    "    f_min,\n",
    "    f_max,\n",
    "    c_approx\n",
    ")\n",
    "\n",
    "# construct circuit factory for payoff function\n",
    "butterfly = UnivariateProblem(\n",
    "    uncertainty_model,\n",
    "    butterfly_objective\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY8AAAE+CAYAAAB1DJw3AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3dd5hU9fXH8fcBRJogKqJGZe01iYpRSTQKFuwELFgRLFhoSzTVXyJqTNEki6gIKIpEEBuKHVHBxC5YYkTsgAiKUqQsKuX8/jh3dRhmd2Z2Z+Y7d+a8nmeeZe7cO/PZYfaeufd+i6gqzjnnXDYahQ7gnHMufrx4OOecy5oXD+ecc1nz4uGccy5rXjycc85lzYuHc865rHnxcEVLRIaIiCbc5ovI/SKyU8BMPxOR10TkaxHRaFlLEZkgIouinL1r2XZM0u9Tc7u1oL/E93n6isgvUiyfLSJ/D5HJxUeT0AGcS+Mr4Ojo3zsCVwNPi8heqroyQJ6RwEKgK/BNtOxi4ASgF/Ap8GEd288C+iQtW5jjjJnqC/wPeDBpeXdgUeHjuDjx4uGK3RpVfSn690siMhf4D3AscG+APLsDo1T12aRl76rq/RlsvzLh9ylKqvp66Ayu+PlpKxc3M6KfFQAi0klEHhKRBSKyUkTeEJEza1YWkc2iU0y9E59EzEciUpWwrIuIvByt/7mIDBeRVtFjh0WnqRoD10enm8aIyGzgPGDfmtNQ9f3Fal5DRPZOWj5NRO5LuD9GRKaLyJEi8t/o935ORPZK2q6xiPxORN4TkW9EZJ6IjKl5TqAjcE7C6bPe0WMbnLYSkVNF5K3oeT4RkWtEpEnC472j5/ihiEyJMs0SkR71fT9ccfPi4eKmIvr5WfSzA/A8tgM/AbgfuF1ETgdQ1cXAA0DvpOc5DNgBuA0g2vE+AXwJnARcAZwB1Oy0XwM6Rf/+R/Tvq7FTPI9hp6M6JaxTKxFpknjL6Lfe0PbAdcA1wOnAlsDdIiIJ64wErgTuAY4HLgVaRI9dEmV+LCH3o7XkPQq4G3sPugE3AJcBN6ZYfTzwEPa+vA9MEJFt6/k7uiLmp61c0UvYwe4IDAeWA08BqOqEhPUE+DewLXABcFf00GjgSRHZUVU/ipb1AWao6lvR/T8Ac4ATVXVt9HyLsR1yJ1V9ETttBjA78dSTiHwBtM/wdFRHYHXS77eLqn6QwbaJNgN+pqrvR8/RCCuSuwGzRGR3rKAOUtVhCdvdDaCqM0VkJfBFBrmvAqap6jnR/Sei9+EvIvInVZ2XsG6VqtYU5BnA51jhGpHl7+eKnB95uGK3ObazXQ28ixWQnqq6AEBE2orIMBGZk7BeX2DXhOd4GisM50TbbIIdXdyesM4BwAM1hSNyP7AGODiHv887wE+Sbp/U43lm1xSOyMzoZ823/M7RzzH1eO7viEhjYD82vL50N7b/SD7SerLmH6q6CGsM4EceJciPPFyx+wo4AlDsVNV8XX8o6DHAQdgppJnAMqz1U7eaFVRVReR24FwRGQKcil27GJ/wPFtj35JJ2G6tiCzCvuXnSrWqTs/B8yxNuv9t9LNZ9HNz7OL8sga+zhbARiS9Nwn3k9+bVLma4UqOFw9X7NbUtrMVkWbYKZF+qjoiYXmqI+rbsesYnbHrHw+q6pKExxdg1w0Sn78xthNe3JBfIAtfRz+bJi1vi12LycYioKWItG5gAfkSO5rbMml5++hnod4bV2T8tJWLs42xz3BNf4uaU1InJq+oqp9gp1SuxE5D3Z60ystA96hg1OiBfcF6Lrexa1Vz7WCPmgUish3WFDhbz0Q/e9WxTtqjgug03gzglKSHTgXWAS/WI5srAX7k4WJLVb8SkVeBP4rIMmxn9lvsVFfrFJuMxs7dzwOmJD32J+B14EERuRk7T/83YHJ0sTzvVHWeiEwHrhaRaqww/p56fLtX1XdFZBTwDxHZEmtIsClwsqqeFq02C+gqIl2xI5WPo+sUya4AJken/iYAP8ROE96SdLHclRE/8nBxdwbwETAWuB67yD22lnUfwS6A36Gq6xIfUNW3gWOw0zMTsWJyF3ByfmLX6nRgLnAn8GespdO79XyuS7AjrbOwJrlDgeqEx/+EXcC/B3gVa+q8AVV9EjgN2B94GKjEmiv3r2cuVwLEp6F15UJEjsUKyK71aBrrnEvgxcOVPBHZBtgF69w2V1WPDxzJudjz01auHPTF+np8DQwInMW5kuBHHs4557LmRx7OOeeyVhZNdbfYYgutqKio17YrV66kZcuWuQ2UR3HKG6esEK+8ccoK8cobp6zQsLwzZsz4UlXbpXxQVUv+1rFjR62vqVOn1nvbEOKUN05ZVeOVN05ZVeOVN05ZVRuWF5iutexX/bSVc865rHnxcM45lzUvHs4557LmxcM551zWvHg455zLmhcP59IZNw4qKji0SxeoqLD7zpW5sujn4Vy9jRsHfftCdTUCMGeO3Qc488yQyZwLyo88nKvL5ZdDdfX6y6qrbblzZcyLh3N1mTs3u+XOlQkvHs7VZfvts1vuXJnw4uFcXS67bMNlLVrANdcUPotzRcSLh3N1mT/ffm6zDQogAjfc4BfLXdnz4uFcbaqrYeRI6N4dPv2U12+6CVQ3vIDuXBny4uFcbcaOhcWLYfBgAJbtuSccdBBcfz2sWxc4nHNhefFwLpV166xIdOwIBx/8/fLKSvjgA3j00XDZnCsCXjycS2XyZJg1y4qFyPfLTzoJttsOhg4Nl825IuDFw7lUqqpg663h1FPXX96kCfTvD888A2++GSabc0XAi4dzyf73P5gyxYpE06YbPn7BBdZc9/rrC5/NuSLhxcO5ZEOHQrNmcOGFqR9v2xZ697Zxrz7/vKDRnCsWXjycS/TFF3DnndCrF2y+ee3rDRwI334LI0YULptzRcSLh3OJRoyAb76xC+V12W03OO44GD7c1neuzHjxcK7GN9/ATTfB0UfDHnukX7+yEhYuhLvuyn8254qMFw/nakyYYNcwok6BaR1+OOy9t10jUc1vNueKjBcP58B2/kOHwp57wpFHZraNiB19vPkmTJuW13jOFRsvHs4BPPssvPHGhp0C0znjDNhiC+806MqOFw/nwDoFbrEFnHVWdts1bw4XXwwPP2zDljhXJrx4OPfBB7bzv+giKwbZuvhi63k+bFjuszlXpLx4ODdsmO38L7mkfttvvTWcdhrcdhssXZrbbM4VqYIXDxHZU0SeFpFqEZkvIleJSOMstm8kItNFREXk+HxmdWVg6VLb6Z92mhWB+qqshJUrYfTo3GVzrogVtHiISFvgKUCBbsBVwKXAlVk8zfnAtrlP58rSrbfaTj/T5rm12W8/+PnPbZbBNWtyk825IlboI4+LgOZAD1WdoqojsMLxSxFpnW7jqPhcA1ye35iuLKxZYzv7Qw+Fffdt+PMNHgxz5sCkSQ1/LueKXKGLxzHAZFVdlrBsAlZQDs1g+6uB54Gn85DNlZsHHoC5cxt+1FHjhBNghx2s5ZZzJa7QxWN3YFbiAlWdC1RHj9VKRH4EnAtclrd0rrxUVcGOO8LxObp01rixDZj4/PPw6qu5eU7nipRoAYdVEJHVwK9UdWjS8nnAWFX9fR3bPgu8rKq/FpEK4GPgBFV9pJb1+wJ9Adq3b99xwoQJ9cq8YsUKWrVqVa9tQ4hT3pBZN5k5k479+vF+//58etJJGW2TSd7GK1fS6dRTWdSpE+/83//lImq9xOlzAPHKG6es0LC8nTt3nqGq+6d8UFULdgNWA5Upls8D/lzHdqcBnwGto/sV2EX34zN53Y4dO2p9TZ06td7bhhCnvEGz9uyp2rq16rJlGW+Scd7KStUmTVTnzatfthyI0+dANV5545RVtWF5gelay3610KetlgBtUixvGz22ARHZCLgO+BvQSEQ2BWourrcUkU3yEdSVsE8+gfvusxkBN8nDx2fgQFi3zkboda5EFbp4zCLp2oaIbAe0IOlaSIKWWNPcf2IFZglQM3n0BOD1vCR1pevGG20gxAED8vP8O+wA3brByJFQXZ2f13AusEIXj8eBrklHCz2BVcCztWyzAuicdDs9euz3wJn5iepK0ooVMGoU9OgBHTrk73UqK2HxYpuV0LkSVOjiMQL4BpgoIkdEF7WHAP/UhOa7IvKBiIwGUNU1qjot8Qa8FK36lqq+XNhfwcXaHXdYr/JcNc+tzSGHWMdBn+vDlaiCFg9VXQIcDjQGHsY6CFYBVySt2iRax7ncWbcOrr8eDjgAOnXK72vVzPXxzjvw5JP5fS3nAmhS6BdU1ZlAlzTrVKR5fDaQxaQLzgGPPQbvv2/TxmYzZ0d99ewJv/619Sfp2jX/r+dcAfmouq58VFXBtttChv06GqxpU+jXDyZPhpkzC/OazhWIFw9XHt58E555Bvr3h402KtzrXnghbLyxz/XhSo4XD1cehg6FFi2gb9/Cvm67dnD22TB2LCxaVNjXdi6PvHi40vf55zB+PPTuDW3bFv71Bw2CVausibBzJcKLhyt9N98M335rPb9D2HtvOPJI65z47bdhMjiXY148XGn7+msYPhyOOw522y1cjspKmD/fhkVxrgR48XClbfx4+OKL/HcKTOfoo614VVV5p0FXErx4uNKlahfKf/hD6FJn16L8a9TIrn1Mnw4vvBA2i3M54MXDla5nnoG33rJTRoXoFJhOr152wX7o0PTrOlfkvHi40lVVBVtuCWecETqJadnSmgpPnAizZ4dO41yDePFwpendd+HRR+Hii6FZs9Bpvtevnx0F3Xhj6CTONYgXD1eahg2z4UEuvjh0kvVttx2cfDLceissXx46jXP15sXDlZ7Fi2HMGDjzTGjfPnSaDQ0eDF99ZRmdiykvHq703HKLzeA3aFDoJKkdeCAcdJAND79uXeg0ztWLFw9XWlavhhtusKa5P/5x6DS1GzwYPvwQHnkkdBLn6sWLhyst990Hn34avlNgOj162PUPb7brYsqLhysdqtY8d5dd4NhjQ6epW5MmNjz81Kk2XLxzMePFw5WOF1+EV1+1ax2NYvDRvuACGybejz5cDMXgL8y5DFVVwaabwjnnhE6SmbZtbZj48eNt2HjnYsSLhysNs2dbz+2+faFVq9BpMjdokA3TfvPNoZM4lxUvHq403Hij9dzu3z90kuzsuqsNF3/zzTZ8vHMx4cXDxd/y5da34+STrQVT3FRWwsKFcNddoZM4lzEvHi7+br8dli0r/ua5tTn8cJttcOhQn+vDxYYXDxdva9daT+1OnazndhyJ2NHHf/8L06aFTuNcRrx4uHh75BH46KP4HnXUOPNM2GILazHmXAx48XDxVlUF228P3buHTtIwzZrZCMCPPALvvx86jXNpefFw8fX66/DsszBggPXYjruLL7bfY9iw0EmcS8uLh4uvqiqbne/880MnyY2tt4bTT7cGAEuXhk7jXJ28eLh4WrAAJkyAc8+1XuWlorISVq6E0aNDJ3GuTl48XDwNHw5r1sDAgaGT5Na++8Khh9qpqzVrQqdxrlYFLx4isqeIPC0i1SIyX0SuEpHGabbZS0SeiNb/RkTmisitIrJ1oXK7IrJqFYwYASecADvvHDpN7lVWwty58OCDoZM4V6uCFg8RaQs8BSjQDbgKuBS4Ms2mbYCPgcuArsAVwBHAYyJSAldKXVbuvBO+/DL+zXNrc8IJsOOO3mzXFbVC73gvApoDPVR1GTBFRFoDQ0Tk2mjZBlT1BeCFhEXTRGQe8CTwI+C1POd2xULVemLvs4+d3ilFjRvb6bjKSnjlFTjggNCJnNtAoU9bHQNMTioSE7CCku2eYFH0s2kugrmYmDIFZs60ow6R0Gnyp08f2GQT6z3vXBEqdPHYHZiVuEBV5wLV0WN1EpFGItJURHYD/gq8CrySj6CuSFVVQfv20LNn6CT51bq1NUG+5x6bVte5IiNawIHYRGQ18CtVHZq0fB4wVlV/n2b7J7BrHgAzgGNVdWEt6/YF+gK0b9++44QJE+qVecWKFbSK0fwQccqbbdYWc+ZwQO/efNynD3N69cpjstQK/d42W7CAA886i7mnncbHF1yQ1bZx+hxAvPLGKSs0LG/nzp1nqOr+KR9U1YLdgNVAZYrl84A/Z7D9LsCBwFnYEcwMoFm67Tp27Kj1NXXq1HpvG0Kc8madtW9f1Y03Vl24MC950gny3nbvrrrZZqorV2a1WZw+B6rxyhunrKoNywtM11r2q4U+bbUEazmVrG30WJ1U9X1VfVlV78SOQPYFzshtRFeUFi2CsWPh7LOhXbvQaQpn8GBYvBj+9a/QSZxbT6GLxyySrm2IyHZAC5KuhaSjqnOAxcCOOUvnitfIkTbT3qBBoZMU1sEHw377WQuzdetCp3HuO4UuHo8DXUVkk4RlPYFVwLPZPFF00XxzrP+HK2XffmvTzB55pE2aVE5E7Ohj1ix48snQaZz7TqGLxwjgG2CiiBwRXdQeAvxTE5rvisgHIjI64f7fReSvItJdRDqLyCXAZOBDrKmvK2X33GNjWZVqp8B0Tj0VttrKjj6cKxIFLR6qugQ4HGgMPIz1LK/CeownahKtU2M6cAgwGngUGAjcDxykqivzHNuFpGrNc3ffHbp2Tb9+KWraFPr1g8mTrY+Lc0Wg4EN7qOpMoEuadSqS7k/AjzDK03PPwWuvwc03Q6MyHsfzwgvhmmus0+DIkaHTOOej6roiV1UFm20GAfp1FJV27eCss6zF2aJF6dd3Ls+8eLji9dFHNrLshRdCixah04RXWWktzvzIwxWBtMVDRHqJyOaFCOPcem64wQYJ7NcvdJLisNde1uLsppusBZpzAWVy5HE7sBOAiKwVER/i0+XfsmU2m17PnvCDH4ROUzwqK2H+fLj33tBJXJnLpHgsAbaJ/i3YXBzO5dfo0bB8ue0s3feOPhp2282a7RZwXDrnkmXS2uop4F8i8i5WOMaISK3NY1XVj0xcw6xda9OwHnww7J96TLay1aiR9bK/5BJ4/nl7j5wLIJMjj3Ox/hivY0ceHwNv13FzrmEmTYLZs8u3U2A6vXpB27beadAFlfbIQ1Wrgb8DiMgRwOWq+ma+g7kyVlUFFRXQrVvoJMWpZUvo2xeuu86KbEVF6ESuDGXS2mqtiPwkujsNSDlVrHM5MX26dQwcONBaWrnU+vWzca9uuCF0ElemMjlt9S2wcfTvXkAZjYftCq6qyqZfPe+80EmK23bbwSmnwK23WsMC5woskwvmM4EhIvIgds3jZBGp7SqmqurNOUvnysunn9ogiP372zSsrm6VlTBhAtx+ux2pOVdAmRSPAcBIbABDBS6rY10FvHi4+rnpJpuzwneEmTnwQOjUyVqm9evnp/lcQaU9baWqL6jqD1V1I+zI4yBVbVTLzT+9rn6qq23YjW7dYIcdQqeJj8pK+PBDePTR0Elcmcl2bKvO2Gks53Jr7FibbtWb52anRw+7/lFVFTqJKzNZDcmuqs8CiMiBwMHAZthUsM+p6su5j+fKwrp11mehY0fv9JatJk1gwAD49a/hjTdgn31CJ3JlIqsjDxFpKSKPAS8Af8E6EP4FeEFEHhURH/rUZW/yZHj3XTvqEAmdJn7OP99GHfZOg66Asj1tdS3QCTgNaKaqWwPNovudgL/lNp4rC1VVsM021vTUZa9tW+jTB+66Cz77LHQaVyayLR4nAb9R1XtVdR2Aqq5T1XuB3wL+1++y87//wZQp1lqoadPQaeJr4EAbpn3EiNBJXJnItni0AT6p5bFPAG+c77IzdCg0b24TPrn623VXOO44GD7cJoxyLs+yLR5vAheLrH9iOrp/cfS4cxnZaOlSuPNOG+hvc59vrMEGD4YvvrDTV87lWVatrYDfA48Ds0TkAeBzYEugO1ABHJPTdK40jRsHl1/OT+fMsfs77RQ2T6no0gW23RYuvJBD16yB7beHa66BM88MncyVoGyb6j4jIvsBf8Cub2wNLABeBnqoqvcBcXUbN85GhK2u5rvD1yFD7IK57+QaZvx4WLgQVq+293bOHHuvwd9bl3PZnrZCVd9W1dNUdSdVbRH9PMMLh8vI5Zdbb/JE1dW23DXM5ZdvOLe5v7cuT7Lt5/EPEdkzX2FcGZg7N7vlLnP+3roCyvbIozvwloi8IiIXiUibfIRyJWz77bNb7jLn760roKyKh6ruCBwBzMJmF1wgIuOjGQadS++aazYc/bVFC1vuGuaaa+y9TNS8ub+3Li/qc81jqqr2ArbChmvfFpgsInNE5EoR2THXIV0JOeggWLsWWrdGRaBDBxg1yi/o5sKZZ9p72aEDWrPsxBP9vXV5kXXxqKGqK1R1NHAF8DywHfA74D0RmSQiHXKU0ZWSYcNgo43gnXd49plnbA5u37nlzplnwuzZPDt1Khx2GDz/PKxeHTqVK0H1Kh4iUiEiV4jIR8CTwAqs6e4mwIlYn48JuQrpSsRXX8Ftt0HPntY01+XX4MEwbx5MnBg6iStB2ba26iUizwAfAOcAtwM7qOqxqnq/qn6jqo8BA4Hapqp15Wr0aFixwiYwcvl3/PGw884+14fLi2yPPEYCnwFdVXVHVb1aVeelWO894E+pnkBE9hSRp0WkWkTmi8hVIlLnDIQi8hMRuV1EPoi2ezc68mmWZX4Xypo1dsrqkENs3g6Xf40awaBB8PLL8OKLodO4EpNt8dgm6hD4dF0rqeoCVb0yebmItAWewuY67wZcBVwKbLBukp7ATtiQ78cCNwG/BMZlmd+FMmmS9Xj2mQILq3dv2HRTP/pwOZft8CRLGvh6FwHNsaFMlgFTRKQ1MEREro2WpfJXVf0y4f40EfkaGCkiHVR1TgNzuXwbOtTmJj/xxNBJykurVnDBBfCPf1jx7uDtWFxuZH3BXER6ishTIjJXRBYm39JsfgwwOalITMAKyqG1bZRUOGq8Hv30K6/Fbvp0eO45m3MiuY+Hy7/+/W2GxhtuCJ3ElZBsL5ifAdyBXTDfFngIeCR6nmXAjWmeYnesg+F3VHUuUB09lo1OwDrgwyy3c4U2dChssgmce27oJOVp++3hpJPg1lth+fLQaVyJyPbI41fA1UC/6P5wVT0X2AH4EisCdWkLLE2xfEn0WEZEZCvg/4B/qWq6ox0X0qefwt13w3nnQWufKyyYwYOtqfSYMaGTuBIhqpp+rZqVRVYAx6vqNBFZDRypqtOix7oDVapaUcf2q4FfqerQpOXzgLGq+vsMMjTFLrpvC3Ss7TqMiPQF+gK0b9++44QJ9et2smLFClq1alWvbUMotrw73Hor248fz8t33snXSX07ii1rOnHKmyrrvv360XTpUl4eO7boTh/G/b0tZg3J27lz5xmqmrrbhapmfAPmY810AWYDFyc81gNYnmb7hcAVKZavxIpKutcX7BrJImD3THN37NhR62vq1Kn13jaEosq7cqXqZpupdu+e8uGiypqBOOVNmfXuu1VB9cEHC54nndi/t0WsIXmB6VrLfjXb01avAj+K/v0Q8EcRuUBEzgGuA15Ks/0skq5tiMh2QAuSroXUYijWxLebqmayvgvpzjth8WLvFFgsevSw6x/ebNflQLbF4y9AzeQAfwReAW7Gepp/CVyYZvvHga4isknCsp7AKuDZujYUkd8B/YGzVPW5LHO7QlO1C+X77WcdA114TZrAgAHw7LPw+uvp13euDhkVDxFpLiInAT8DmohIe1VdqqrdgJbApqp6oKp+lOapRgDfABNF5IjousQQ4J+a0Hw36kk+OuH+GcCfgbHApyJyUMKtXRa/ryuUJ5+Ed96xow6R9Ou7wjj/fGjZ0gq7cw2QtnhEQ6y/DdyLnZr6F/CuiBwFoDaeVW2d+9ajdnH7cKAx8DDWs7wKG5k3UZNonRpHRT97Ay8m3Y7L5LVdgVVVwVZb2SCIrnhsuin06QN33QULFoRO42IskyOPa7H+FIdg1yb2wjrojazPC6rqTFXtoqrNVXVrVf2Dqq5NWqdCVXsn3O+tqlLLbUx9crg8mjkTJk+Gfv2gadPQaVyyQYNsrLHhw0MncTGWSfHoBPyfqj6vql+r6jvYtY3tRWTr/MZzsTRsGGy8MVyY7hKYC2LnneGEE2DECFi1KnQaF1OZFI+tgeRrGR9izWa3ynkiF2+LFsHYsXD22dDOL0cVrcGD4csvrUWcc/WQaWurzHsSuvI2apR9mx00KHQSV5dDD4V99rEL51l0FHauRqbFY3LS4Ic1V9qeznJgRFfKvv0WbrwRjjwS9t47dBpXFxFrCTdzJkyZEjqNi6FMhmRPN9eGc+a++2D+fLjlltBJXCZOOw1+8xtrGXfUUenXdy5B2uKhKSZ1cm4DqrYT2m03OPro0GlcJjbe2FrE/fGP1idnjz1CJ3IxkvV8Hs6l9MILNm/HoEE2/amLh4susiLinQZdlvyv3OVGVRW0bQu9eoVO4rLRrp21jBs71lpfOZchLx6u4WbPhgcegL59begLFy+VlfD119ZSzrkMefFwDXfDDdZ6p3//0Elcfey1l7WQu/FGazHnXAa8eLiGWb7cpjc95RTYdtvQaVx9DR5sY13dc0/oJC4mvHi4hrn9dli2zHY+Lr66doXdd7drV95p0GXAi4erv7VrbRyrTp3ggANCp3EN0aiRXft47TX4z39Cp3Ex4MXD1d8jj8CHH/pRR6k4+2zYbDOfadBlxIuHq7+hQ21a0+7dQydxudCihfX7mDQJPko3r5srd148XP288QZMm2bTmjbJZJQbFwuXXAKNG9vpSOfq4MXD1c/Qodan47zzQidxufSDH9jsj6NHw1dfhU7jipgXD5e9zz6zaUx797Ze5a60DB4MK1ZYAXGuFl48XPZuvtk6k/mcHaWpY0c45BA7dbVmTeg0rkh58XDZ+fprKx7HHw+77BI6jcuXwYNhzhx48MHQSVyR8uLhsjN+PHzxhTfPLXUnngg77OCj7bpaefFwmVO1ncmPfgSdO4dO4/KpcWMYOBCefx5efTV0GleEvHi4zD3zDLz1lvVEFgmdxuXbuefCJpt4p0GXkhcPl7mhQ2HLLeH000MncYXQujWcfz7cey/Mmxc6jSsyXjxcZt57z4YjufhiaNYsdBpXKAMHwrp1Nly7cwm8eHgf2NUAAB9lSURBVLjMDBsGTZta8XDlo6LChp8ZNQpWrgydxhURLx4uvSVLbOj1M86A9u1Dp3GFNniwfQbGjg2dxBURLx4uvVtvhepqu1Duys9Pfwr772/XvNatC53GFQkvHq5ua9bYNLOdO8OPfxw6jQtBxI4+3nsPHn88dBpXJLx4uLpNnAiffOJHHeXulFNs0ERvtusiXjxc3aqqYKedbDgSV7422gj694enn4b//jd0GlcECl48RGRPEXlaRKpFZL6IXCUijdNs01RErhOR/4jIKhHxSZYL4aWX7DZokE1T6spb3742YdT114dO4opAQfcIItIWeApQoBtwFXApcGWaTVsA5wPVwAv5zOgSXH89tGkDffqETuKKwWabwTnnwLhxsHBh6DQusEJ/nbwIaA70UNUpqjoCKxy/FJHWtW2kqkuBzVS1K/BAYaKWuU8+sZ7F558PrVqFTuOKxaBB8M03NrKyK2uFLh7HAJNVdVnCsglYQTm0rg1V1U9VFdJNN9lAiAMGhE7iisluu8Gxx8Lw4TY8vytbhS4euwOzEheo6lzsdNTuBc7iarNypfUo7tEDOnQIncYVm8GD7bTVXXeFTuICalLg12sLLE2xfEn0WM6ISF+gL0D79u2ZNm1avZ5nxYoV9d42hFzk3WbSJHZdsoTXDj2UZXn83cvxvS2UvGZt3Jj9d9wR/vQnpldU5GSEZX9v8ydveVW1YDdgNVCZYvk84M8ZPkd/orNYmd46duyo9TV16tR6bxtCg/OuXau6666qP/mJ6rp1OclUm7J7bwso71lHj1YF1aefzsnT+XubPw3JC0zXWvarhT5ttQRok2J52+gxF9oTT1hPYp+zw9XljDOgXTvvNFjGCl08ZpF0bUNEtsOa4s5KuYUrrKoq2GYb61HsXG2aNbMRlh95xL5suLJT6OLxONBVRDZJWNYTWAU8W+AsLtn//gdPPWU9iTfaKHQaV+wuucSG6fdOg2Wp0MVjBPANMFFEjoguag8B/qkJzXdF5AMRGZ24oYgcIyInA/tE90+Obt4cKFeGDoXmza0nsXPptG9vp6/GjIHFi0OncQVW0OKhqkuAw4HGwMNYB8Eq4IqkVZtE6yS6GbgXOC+6f29065yvvGXliy/gzjuhVy/YfPPQaVxcVFbacP233BI6iSuwQjfVRVVnAl3SrFORyTKXQyNGWM/hQYNCJ3Fx8uMfQ5cuNk3tL3/ppzvLiI9256xoDB8ORx8Ne+wROo2Lm8GDYd48uP/+0ElcAXnxcHD33fDZZ7YTcC5bxx4Lu+xiLfV8FKGy4cWj3KnahfI994QjjwydxsVRo0Z2uvOVV+DFF0OncQXixaPc/fvf8Prr3inQNcw558Cmm3qnwTLixaPcDR1qravOOit0EhdnrVpZE++JE2H27NBpXAF48ShnH30EkybBRRdZ/w7nGqJ/fzt6vfHG0ElcAXjxKGfDhkGTJtZT2LmG2m47G9bmlltg+fLQaVyeefEoV8uWwW23Qc+eNpaVc7lQWWmfrdtvD53E5ZkXj3I1erR9O6ysDJ3ElZIDD4ROnWy8q7VrQ6dxeeTFoxytXWunrA45BDp2DJ3GlZrBg+162sMPh07i8siLRzmaNMlaxPhRh8uH7t1t+mJvtlvSvHiUo6FDoaICunULncSVoiZNYMAA60P02muh07g88eJRbmbMgP/8BwYOhMbJAxc7lyPnn299P4YODZ3E5YkXj3IzdKj9UZ97bugkrpS1aQN9+sCECbBgQeg0Lg+8eJST+fPtj/m88+yP27l8GjQI1qyBm24KncTlgRePcjJ8uLW0GjAgdBJXDnbaCU480eaKWbUqdBqXY148ysWqVfZH3K2b/VE7VwiDB8OiRTZLpSspXjzKxZ132h+xN891hfTzn8O++9q1Np/ro6R48SgHNXN27Luv/TE7VygidvQxcyY8+WToNC6HvHiUgylT7I/X5+xwIfTsCVtt5Z0GS4wXj3JQVWV/vD17hk7iylHTptCvH0yebF9iXEnw4lHq3nkHnnjChl3feOPQaVy5uugiaNbMOw2WEC8epW7YMCsaF10UOokrZ1tsAWefDf/6F3z5Zeg0Lge8eJSyRYvgjjtsitl27UKnceWushK+/hpGjgydxOWAF49Sdsst1r/Dm+e6YrDnntC1q/U4//bb0GlcA3nxKFWrV9tc0kccAXvvHTqNc6ay0sa6uvvu0ElcA3nxKFX33Qeffmpt7J0rFl27wh57WAtA7zQYa148SpGq/XHuthscfXToNM59T8SOPl5/3eb7cLHlxaMUvfgivPqqjWrayP+LXZE5+2zYfHNvthtzvmcpRVVV0LYt9OoVOolzG2re3JqOT5oEH34YOo2rp4IXDxHZU0SeFpFqEZkvIleJSNop7USkjYjcLiJLROQrERknIpsXInOcbPzZZzBxIvTtCy1bho7jXGqXXGLT1Q4bFjqJq6eCFg8RaQs8BSjQDbgKuBS4MoPN7wEOA84HegM/AR7MR8442/aBB+y8cr9+oaM4V7tttrHhcm67Db76KnQaVw+FPvK4CGgO9FDVKao6AiscvxSR1rVtJCKdgKOAc1T1flV9ADgLOFhEjshL0nHjoKKCQ7t0gYoKu1/Mxo2D7bdn23vusR7lfjHSFbvBg2HFivj8ncVxn5DHvIUuHscAk1V1WcKyCVhBOTTNdp+r6nd7RFV9Bfg4eiy3xo2z0z5z5iCqMGeO3S/WD0tN3k8+QQCqq4s7r3Ng4641agRLlxb/31lc9wl5zCtawLbWIrIQGK6qQ5KWrwSGqOp1tWx3D7Clqh6WtPxRAFU9rq7X3X///XX69OmZB62osDc7WZMmsOuumT9Pobz3ns0VnaxDB5g9u+BxMjVt2jQOO+yw0DEyFqe8scgap7+z2v7GijEr5GyfICIzVHX/VI81qW+2emoLLE2xfEn0WH222zHVBiLSF+gL0L59e6ZNm5ZxyEPnziXVrBe6Zg1fFOEYUe1mzkydd+5cns3i9y60FStWZPX/Elqc8sYha5z+zmr9GyvCrFCgfYKqFuwGrAYqUyyfB/y5ju2mAA+mWH4n8EK61+3YsaNmpUMHVetqt/6tQ4fsnqdQ4pY3MnXq1NARshKnvLHIGqfPbZyyquYsLzBda9mvFvqaxxKgTYrlbaPHcr1d/VxzDbRosf6yFi1seTGKW17nIF6f2zhlhYLkLXTxmAXsnrhARLYDWkSPZbxdZPc029XPmWfCqFHQoQMqYucJR42y5cUobnmdg3h9buOUFQqTt7ZDknzcgN8Bi4FNEpZdBlQDrevYrhPWN+TghGX7R8uOSPe6WZ+2ShCLw/8Eccobp6yq8cobp6yq8cobp6yqDctLEZ22GgF8A0wUkSOii9pDgH9qQvNdEflAREbX3FfVF4EngbEi0kNEfgGMA55T1acK+hs455wrbPFQ1SXA4UBj4GGsg2AVcEXSqk2idRL1BJ4FbgPGAjOA7vnM65xzLrVCN9VFVWcCXdKsU5Fi2VKgT3RzzjkXkI+q65xzLmtePJxzzmWtoMOThCIiXwApxkHIyBbAlzmMk29xyhunrBCvvHHKCvHKG6es0LC8HVQ1ZRf6sigeDSEi07WWsV2KUZzyxikrxCtvnLJCvPLGKSvkL6+ftnLOOZc1Lx7OOeey5sUjvVGhA2QpTnnjlBXilTdOWSFeeeOUFfKU1695OOecy5ofeTjnnMuaFw/nnHNZ8+LhnHMua148nHPOZc2Lh3POuawVfFRdlxvRDIzHAgLcq6qLRGRbbHKtnYDZwChVfStcShCR3wCPhc6RKRFpDjRR1eUJy9oB/YE9gXXAG8BwVf0qTErnwvOmuhEREWx+kOOAPYDNsB3FZ8BLwBhVfS9cwu+JyAHAFKAlsAabnbEr8BiwFngb2BvYCptp8T+BoiIi67AZH2cB44EJqvphqDzpiMhjwPuqOii63wl4HPsszMCKdUfgW6CLqr4dMOu+QHNVfSFh2dHYjJ01he5NYEjiOsUi+ps7AdgP+4xMx75oFPVOSURaY2NFdVHV50Lnge8ydQGaAo+q6sroS08/YEfgI+zL5PycvWaR/z8VRPQmP4btFD7Ddgw/wD7Qj2Nv/m7A1ap6daicNURkCnbU2B1YiU2o9Qts53ayqq4WkY2BB4Fmqto5YNZ1wN+AHwJHYrlfwwrJPar6aahsqYjIl8B5qjopuv8S9h7/ouZoRETaAA8BX6tq14BZXwIeVtVrovvnArcCU4FnsEJ3OHAIcFLN7xQo6wvY+/pOdL8tNjtoR2BFtFor7Ita18QjvxBE5JI6Hm4OXAdcD7wPoKrDC5ErFRHZGXga2C5a9DFwFPYFc1PgQ2z/tQroqKrzcvLCtc1PW0434C7sQ/DDhGXbAE8A90f3D8U+5OcWQd5FwDEJ97fEvmUelbTeccCXgbOuAw6I/r0Z0Df6oK+JbtOiZZuFfl+jjNXAzxPuf5v8via8tysDZ12WmA34ALghxXojgDeL5XMQ3R+NHTEfnbDsaGAJUFUEn4N12FH8ulpuiY+tDZz1HuwIc+fob+xf0f7sBWCTaJ0tonVG5up1/YK5OQb4rSacl1c7vLsI+IWIbK2qzwJ/BgYFyphMU/w7+TCyqA4rVXWxqo5S1cOBbYFLscPsEcACEXk0aEDzPyDxSO1z7A8y2eZYoQlpXdL9DsB9Kda7D/vmWUxOBK5S1SdqFkT/vgboESzV9yYBC4HzgMaq2qjmhn0eBDgsWpY8ZXahHQxco6ofqOpi4P+w655/1+gITlW/BIay/me7Qbx4mEbYN4lka7EPSZvo/svAroUKVYfpwGUi0kpEGgG/Bz4FLhaRxgAi0gS4BNsZFh1V/UxVr1fVnwI7YPPYbxM4FsBfgd+KyLnRe3gNcJ2IHCkiTUVk4+i6wl+wb3wh/Qc4M+H+20Cqobd/gn0+ismm2Oc42QzsWl1QqtodOAf4FfCqiPws8eEwqWrVFjvdXqPm/zp5DqOPsC9tOeGtrcwU4E8i8l9V/Qi+Oyc7DPtPqblQ3goohhY2l2OZl2Cnfqqxi2X3Au+LSM0F822wUwFFTVXnYDvtvxZBlokiMgD7llYFvIt9eXgC22lItOpD2I4lpN8Dz0dfIG7ALpTfISKbYacDBftcVAK/DRUywUkiUlPclgCpJhnaAjsdF5yqPikiP8Lev0dF5AmsNWPQ6zEpLMSOOmusBUZiR82JtiSH2f2CORA1cZ2MHVXMwc5z7wB8A5yuqo9H612LzazVM1TWGlHm47EvAPer6gIR2Qr4NXaKYg5wq6q+FjAmInIFcIvmsJVHIYjI5kBP4ADsm7BgO7x3gEdUdUbAeN8RkX2Am4EDWb+41fx7CXZ66PowCU3UcCLZGFU9N2m9kcCeqnpIYZJlJvrbuhY7pTYSKyidVfXfQYMBIvIgsDj5vUyx3g3A7qp6ZE5e14uHiU73nAr8GGiGXXwcH51DdK6oicgeWAFJLnQvqOrqkNmyISIXAB+q6jOhs6QSNd2uwr6gHadF0ARaRNoDLVT14zTr/RJrOPF0Tl7Xi0fpEZFGqprqm17REJFm2EW9dcAHxbiDi6557EhCnx9VnRs2lXPFwS+YJxGRvUTkJBE5P7qdJCJ7hc6VTER6iMiDIvKYiJwQLespIrOB1SIyJ/oWF5SInBX1P6i530RE/oo10/wvdkF/sYgUwzl5AESko4g8hJ0ffgd4HngR+FhEPhWRq0SkRdCQJUQioXOkIiLNk/+vRWSfaL/QMVSuohCyfXIx3YBzsesEqdp2r8WG++gTOmeU9dQo13NYk8Jq4ALsWs1orFfpXVHuroGzzgQuTrj/jyjvH4CfYc0Mh2AdmH5fBO/tUdi1rulY0+whWEfRb6PMl2Ktmt4A2hZB3uOxfjPvRJ+Fn6dY50DC90U4iqjPQcKyX2AdRtcAq6P3/LjQ72mUrQ3wQJRrDXAL0Bi4I2m/8BywRei8Gf5OJ+XycxD8FyqGGzAg+pDchPXG3SL6oDSO/n0wcGO0U+lXBHlfBUYk3D8zyvaPpPVuB54KnLUaODTh/kJgUIr1LgPmFMF7OwO4o5bPyGzsaL1ZtNMbHjjrkdEO7Pno8zkjuv8PolPS0XrFUDzWsn4nwe7RDviF6P/+Uuzobg0pOmUGyDsMG4JkANAr+sJwP/BJVAjbYf3DPgVuDp03w98pp8XDr3kAIvIRtjO+Ns16vwYuUtUdC5Os1hzLgB6q+lR0vw12gfQITbjQGJ3OGqmqwfpPiMgCoL+q3h/d/wY7GpqWtN6RwEOq2rzwKdfLsQo4UVWnJC1vi/Xs30tV3xGRXsDfVHXrEDmjTM9h43D1SVh2Lrbjm4K1FPxaRA7ELpwH68wWtbY6SFVfie6/BnyqqickrfcY0FJVDw0QMzHHbKzj3S3R/X2x4txHVe9IWO8C7Ih5hyBBLcNtGa7aAevYmJPPgV/zMFsDr2Sw3isUQQcmrBlm4gegZmygpUnrrcA6Y4X0ENahsWl0/yng9BTrnY59uwttIdbiLtmPsfe9pp/PHL7vPBrK3sCdiQtU9TZsKJ2DgGeiPh/FaG+syWuyUdhAiaG14/v+XRCNYYWNE5XoA1L3Vymkc7CjoR+muXWo7QnqwzsJmjeBC0Tk31pLK6Xogt4F2EXe0OZgo6ZOBlDVtVETwneS1tsJWFDgbMl+h/WE/p+I3Ao8DPxNRPbm+45snYF9sRFWQxsFXC0irbBC9y3WQ/tyYKp+319lRyB0y6uvsZGV16OqM6Ie0ZOx00JDCpyrNomnOb4idYe1lRTHl9qPsSL8bHT/EOw020+x6xw1fkb4z8H7wCuq2quulUTkZODuXL2oFw9zKdaDeKaITMSGD6/5Ft8G2B07R7stxdFjeyJJwwyo6ssp1juN9T/oBaeqi0XkIGzn+0uslytAp+j2LXaK5RBVfTVMyu+p6jXRKZbfAn+sWYw1QKhMWHU1dkE9pP9i590fSn5AVT+KCshjwJgC56rNZBFZE/27DfaF4dmkdXYn/BcesPHWrheRH2KF7lTsi9Afoi8Wb2JHSIOB0CNtv4QVtXQSO5E2mF/ziIjITljv7KP5fmjjGp9gLW6u0yKeiyKZiGwPLFXVohjuAUBEKli/I9uHWpx9PDbCjtyaAR8V03tYQ0QuxIYo2Vdr6cwqIi2xVkNHqA3qF0Q00kCy91V1fNJ606LlxdDMfCB2OnUjbLSGESJyOnZNqWZgzFHAb0J+hqMmwz9T1WFp1tsCu2aXXLDr97pePDYUteuuuVawVFVDj57qnCsS0SnsLVT1i9BZQvLiUWKiQ+rXgDOL4TSQxHBaV4nJFL/OheTFI0G009gGO5XyZYrHtwCOVdWxBQ+3fo5j63i4JXZR7LdEw7Gr6mOFyJWKxGhaV4jXFL+Zisa9OkVVrwqco+BTpeZSdMSROG3uDOz3CL4TFRut+CTs72mMqs4SkR8DV/L9F54bVXVyzl40dMeVYrgBG2PDma+NbquxntptktYL3tkqyhGnWc6+BLol3H8J6xG9ScKyNtiF08lF8N5OwaZx3RQ7130jMA/rvb1Rwuflcaz1VfDPbwa/U047h9Uzw85YK8Gaz+WH2E7tI6xAv4oNxf45sG0RvGcvAHsk3G8bZVwX5VzG950cNwmVM8rWFfvy9Vn0vi7DWjAuifLdFP3drcWmU87N64b+TyqGG9aqZinWFHd/YGD0IX4f2CVhvWIpHtOxFil9sLbbibcfRR/qU2uWBc4am2ldoxxxmuJ3+wxvF4X+3BJoqtQG5I3NtLnYCAP3YjMegjWiWAKMTlrvX8BLOXvd0P9JxXDDmub2T1q2FfBv4AugU7SsWIqHYPN+L8SGTNgh4bE20Qd/gzGOAmV9Bbgi4f4nwGkp1usFfFEEeRcl7SDaRe/nkUnrHVsExaPmKDPdrRiOQOcDpybc7xDl6pG0Xh/gvSL4HCQXjy+AyhTrBR9WB2tKfETC/bZR/i5J6x2FNQDKyet6Pw+zHUmd/1T1MxE5HKvWT4nImRRH+3PUPgmjROQe4E/AW9FEL38KmyylvwLjROQTYCzfT+u6CDtVVdNJsBimdYXvp/h9DjtqSpzi9xm1DpnFMsXvcuAZ4NY06x2MNUMPKchUqTlUzNPmrmL9zqI1/04e6qcF1rE0J7x4mPnALtiRxnfU2m6fJiLXY4eFQS+UJ1PVpUB/ERmFtT1/H/gbRTTHssZrWleI1xS/r2DX5R6ta6Vo7pTQgkyV2kBxmTb3eeCPIvJ+lOXv2GjWv4lGzVgejX/3a6zY5YS3tuK7gcV2VNXD6ljnd9i3ZtWAA8zVRUROw6bK3BYbAC34FJk1JCbTukKspvj9A9BXVZM7tSav93PgSlXtXJhkKTMEmSq1viRG0+aKyM7YUDo1n4PZ2NH8fViP/TlABfZlqLOqvpGT1/Xi8V0zt57AX7SOaWdF5Azs3Hef2tYJLTql0hJYoaprQ+dxDsJNlZpvUiTT5kb9u36GtRB8WlVXRZ2dz+f7LzzjVXVezl7Ti4dzzrlsFcPolS5PROQWERkdOkcm4pQV4pfXuVzzC+ZZEJFbgEaqel7oLBnqTHy+IMQpK8Qor4g8hZ1lODx0lnTilBXilTfXWb14ZCc2OwwAVd05dIZMxSkrxC6vEJ/PbZyyQrzy5jSrX/MoYVETzS1VNfRkNWnFKSvEL69zuRaXilkURKRZNEdGXByHzYgWB3HKCjHKKyIbxeVzG6esEK+8uc7qxSM7sdlhuPIgIv1E5EMRWSUib4rI2SlW248i+NzGKSvEK2+IrH7NI4ZEJNM25al6xBZUnLJCvPJGnUJvwKbIfR2binSMiHQDzlLVnA1F0VBxygrxyhsqqxcP4rXDiPwcG+ZjZpr1imFYijhlhXjlvQz4u6p+N25VNB7bOGCqiByvqouCpVtfnLJCvPIGyeoXzAERWUNmO4wfAAeGHp5ERN4EZqlqzzTrnQzcHTJvnLJGOWKTV0SWAyeo6rSk5RXYfCONsfG32gEveNbMxSlvqKx+zcO8DfxPVU+p6wb8M3TQyEvAQRmslzjwYChxygrxyvsVNjDfelR1Nnbq4kvgReAnhY2VUpyyQrzyBsnqRx58N7jZ0araIc16J2FzWgctuiKyE7CXqj6UZr3mWHPS5GGvCyZOWaMcsckrIpOA5ap6Vi2PN8cGxzuGwAN6xilrlCc2eUNl9eJBvHYYztUQkVOAwcDxtQ3oKSKNgZuxAT13KGS+pByxyRpliU3eUFm9eDjnnMuaX/NwzjmXNS8ezjnnsubFw5UVEektIjNEZLmILBGR10UkL63oRGRXERkiIptmsO4QEdGE23wRuT+6Hpdu297RNq1yk9y59Lx4uLIhNpXwrcBkoAfQC5gEnJinl9wVuAJIWzwiXwGdottlwD7A0yLSMs12j0bbVNczp3NZ8x7mrpz0B0aq6u8Tlj0sIleGCpRkjaq+FP37JRGZC/wHOBa4N3nlqAVNY1X9AviicDGd8yMPV142BT5LXqgJTQ5FpCI6BXSGiPwrOr21UESuSN5ORLqIyMsi8rWIfC4iw2tOHYnIYcDD0aofR885O8u8M6KfFdFzjhGR6SLyCxF5G/gaODDVaSsRaS4i14rIHBH5RkQ+FpG/JOU/X0Tejh6fIyK/xrkM+ZGHKyevAQOib/SPpBnv5zrgEeBkbLyrK0TkS1W9CUBE9gKeAKYAJwHbAX8FdsSGgniNaMwh7BTZAuCbLPNWRD8/S1p2LXBVtPxjYL3rIiIi2Om4TsDVWBH6AXBIwjq/Av4cPdc0oCNwtYhUq+qNWeZ05UhV/ea3srgBPwI+woYWWYcNS3MV0DphnYro8SeTtr0F+BSbhhhgAvA+dtqoZp1To207RfePj+5XZJBtCDaMRJPotiswFVgGbB2tMyZ6vn2Stu0dLW8V3e8a3T+xltdqDawArkhaXlOQGqfL6ze/+WkrVzZU9b/AHtgF8uHY2FR/AKanaKn0QNL9icA2wLbR/QOAB1R1bcI69wNrgIPrGXFzYHV0exc7iumpqgsS1vlUVd9I8zxdgMVa+4gJnYCWwL0i0qTmBjwDtOf739G5WvlpK1dWVPUb7FrEwwAich7WAus84PqEVRcmbVpzf2tgbvTz86TnXisii4DN6hnvK+AI7KjhM2C+qiYPAfH5BlttaHPsNFltagbRe7uWx7cDfAgeVycvHq6sqepoEbkW2D3poS1rub8g4ed660StnzYHUo4vlIE1qjo9zTqZjCe0CCtutanJdzypi9G7GbyGK3N+2sqVDRFJLgiISDugDRvuRLsn3a+56D0vuv8y0D0qGInrNAGei+5/G/0s9MRRTwObicjxtTz+IrAK2EZVp6e4LS9cVBdXfuThyslb0fDVT2KnoTpgLaKqgTuS1t0rGqr/fqy11XnAIFVdFz3+J2zKzwdF5GbsOsHfgMmq+mK0Ts03+AtFZAJQrapv5edXW88UrCPkeBG5Cmv5tTXwc1W9UFWXisgQ4HoR6QD8G/siuSvQWVWTC6dzG/Di4crJVUA3YBh2XeIz4AXsovTHSev+Gjutcz/Wn+Jq4LsmrKr6togcgzV3nYi1iror2q5mnTkichkwEBiAHbVU5OMXS6SqKiLdo8yV2Axy84HxCetcKyLzsaG8L8V+x/eAu/Odz5UGH5LduQTR1J0fY9N6PhI2jXPFy695OOecy5oXD+ecc1nz01bOOeey5kcezjnnsubFwznnXNa8eDjnnMuaFw/nnHNZ8+LhnHMua/8P8jm630Wa/TMAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# plot exact payoff function (evaluated on the grid of the uncertainty model)\n",
    "x = uncertainty_model.values\n",
    "def payoff(x):\n",
    "    if x <= strike_price_1:\n",
    "        return 0\n",
    "    elif x < strike_price_2:\n",
    "        return x - strike_price_1\n",
    "    elif x < strike_price_3:\n",
    "        return 2*strike_price_2 - strike_price_1 - x\n",
    "    else:\n",
    "        return 0\n",
    "y = [payoff(x_) for x_ in x]\n",
    "plt.plot(x, y, 'ro-')\n",
    "plt.grid()\n",
    "plt.title('Payoff Function', size=15)\n",
    "plt.xlabel('Spot Price', size=15)\n",
    "plt.ylabel('Payoff', size=15)\n",
    "plt.xticks(x, size=15, rotation=90)\n",
    "plt.yticks(size=15)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "exact expected value:\t0.2598\n"
     ]
    }
   ],
   "source": [
    "# evaluate exact expected value (normalized to the [0, 1] interval)\n",
    "exact_value = np.dot(uncertainty_model.probabilities, y)\n",
    "print('exact expected value:\\t%.4f' % exact_value)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Evaluate Expected Payoff"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "# set number of evaluation qubits (=log(samples))\n",
    "m = 5\n",
    "\n",
    "# construct amplitude estimation \n",
    "ae = AmplitudeEstimation(m, butterfly)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "# result = ae.run(quantum_instance=BasicAer.get_backend('qasm_simulator'), shots=100)\n",
    "result = ae.run(quantum_instance=BasicAer.get_backend('statevector_simulator'))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Exact value:    \t0.2598\n",
      "Estimated value:\t0.2290\n",
      "Probability:    \t0.7939\n"
     ]
    }
   ],
   "source": [
    "print('Exact value:    \\t%.4f' % exact_value)\n",
    "print('Estimated value:\\t%.4f' % result['estimation'])\n",
    "print('Probability:    \\t%.4f' % result['max_probability'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZgAAAEPCAYAAAB/WNKuAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAeAUlEQVR4nO3dfbQdVZ3m8e/DixAaiFEgINJEEE3DQL9wRTKNciPvoUcwgGHQ6RUbjToqdK9Ig4gQ0GEJDogudAFLOwyjJtJCM40QQoDcQADFIEGUBAwakBdBui+JIRAJ+c0fu0IqlXPPqfNS5+bcPJ+1zjqndu3aZ9dO5fxu1a69SxGBmZlZp2013BUwM7ORyQHGzMwq4QBjZmaVcIAxM7NKOMCYmVklHGDMzKwSDjBmdUiaIWkg+zwgaUaT2/dLimJZQ+S9WdIjddZfKeklSduV/O53SgpJxzZTZ7NOcYAx23zMAv6LpP2LKyRtDZwM3BgRa7peM7MWOMCYbT7+H7Aa+O811k0ExpKCkFlPcIAxa5GkCZL+XdJzkl6WtFjSR1otLyJeBm4GptRYfSrwAnBX9t17Spop6beSXpH0uKQLJW1bp77bZJfMPlVI/4qk3xfS9pb0Q0mDklZLmiNpv1b3zbZM2wx3Bcw2ZxExI/e5v7B6b+Be4CrgVeBvgZmS1kXErGybAUDFsuqYBUyRdHBEPAiQBY3JwPcj4vUs367Ai8A/Ai8B44ELgF2AzzS5mxuRtEu2X88D07J9OxeYJ+ndvkRnZTnAmLUoImav/yxJwN3A24FP0PqlrDmkgHEq8GCWdgwwJl9mRCwGFue+/17gFeAqSWdGxNoWvx9gOrAdcEREvJSVfx+wHJgKXN1G2bYF8SUysxZJGiPpm5KeBF7LXtOAd7VaZkT8CbgR+HAWtCBdMnsSuD/33VtJmi5piaRXsu/+P8AoUpBrx5HAXGBVdlltG2AF8HOgr82ybQviAGPWumtJP/5fA44G3gP8C7B9m+XOAv4cmCBpe+AEYHZsPPX5dOAS4F+BDwKHAGdk69r9/l2Aj7AhaK5/vR/Yq82ybQviS2RmLch++P8O+ExEXJVL78QfbfNJ/R+nAnsAO7HpJbdTSEHn/Nx3H9Sg3NeBtcCbCuljCsv/CTwEXFyjjJUNvsPsDQ4wZq3ZjnQF4I0Ob0k7kc4m2nrIUkS8Lul6UhDZE1gSEQ8Xso3Kf3em7h1sERGSngH+IlfnrYEjClnvJJ01PeIOfWuHA4xZCyJihaSfAedLWgmsA84h9VXs3IGvmAV8DvgQ6e6wonnApyUtAn4D/D0wrkS5/wZMk/QwqV/nE8AOhTz/GzgNuEvSlcCzwO7A4cBARFzf9N7YFsl9MGatO430434d8A3ghuxz2yLiftJdW6L2HWkXANeTLmPNAl4G/qlE0eeTbiK4GJgJLKJQ54h4ATgUWAZcAdxO6u/ZCRhyKhuzInX7kcmS3gmcBUwADgDuqTG+oNZ2o0kH+4mkwPhj4IyI+I9CvhOArwD7kf7zXxgRP+zkPpiZWWPDcQZzADAJeAx4vIntrgf6gY+T7sV/D3BTPoOkw0h/Rc4HjgNuAWZJOrrdSpuZWXOG4wxmq4hYl33+EbBLozMYSROA+4DDI+LuLO0Q4KfAURFxR5Y2F9g2Ij6Q2/ZWYOeIOKyK/TEzs9q6fgazPrg06Tjg+fXBJSvnAeC32TqyKcwnks508maTxhOMbq3GZmbWil7p5B8PLK2RviRbB7AvsG2NfEtI+9ny6GozM2ter9ymPIY0P1PRILBPLg818g0W1m9E0jTS9B6MGjXq4L326sxA5XXr1rHVVr0Sv4eP26kct1M5bqdyOtlOjz/++IsRsWutdb0SYCoTEdcA1wD09fXFokWLOlLuwMAA/f39HSlrJHM7leN2KsftVE4n2ymbi6+mXgn1g0CtPpQxbDhDWf9ezDemsN7MzLqgVwLMUjb0teTl+2aeIE3IV8w3njTKuplbos3MrE29EmDmALtn41wAkNRH6n+ZA5DNmTSfNH9T3hTg/ohY0aW6mpkZw9AHI2kH0kBLSBP57Szp5Gz51ohYLWkZsCAiToc0bYak24HrJH2edEZyCbBw/RiYzJeBAUlXkAZhTspex1a+Y2ZmtpHh6OTfjfQMi7z1y+8gzb+0DbB1Ic8U4Ouk5228MVVMPkNELMyC1VeAT5PGyZwWEbd3sP5mZlZC1wNMRCwne0Z5nTzjaqS9BHwse9Xb9iYKU8iYmVn39UofjJmZ9RgHGDMzq4QDjJmZVcIBxszMKuEAY2ZmlXCAMTOzSjjAmJlZJRxgzMysEg4wZmZWCQcYMzOrhAOMmZlVwgHGzMwq4QBjZmaVcIAxM7NKOMCYmVklHGDMzKwSDjBmZlYJBxgzM6uEA4yZmVXCAcbMzCrhAGNmZpVwgDEzs0o4wJiZWSUcYMzMrBIOMGZmVgkHGDMzq4QDjJmZVcIBxszMKuEAY2ZmlXCAMTOzSjjAmJlZJRxgzMysEg4wZmZWCQcYMzOrhAOMmZlVwgHGzMwq4QBjZmaV2Ga4K2Bmmxp3zi0bLU8/cC1Tc2nLv3p8t6tk1jSfwZiZWSUcYMzMrBIOMGZmVgkHGDMzq0TXA4yk/SXdKWm1pGclXSRp6wbbzJAUQ7y+kMt37RB5xle/Z2ZmltfVu8gkjQHuAB4FTgD2BS4jBbrz6mz6HeC2QtqJwNnAnEL6UuBjhbTlrdXYzMxa1e3blD8FjAImR8RKYJ6knYEZki7N0jYREU8DT+fTJH0JWBoRiwvZX46In1RQdzMza0K3L5EdB8wtBJLZpKBzeNlCJL0VOAqY1dnqmZlZp3Q7wIwnXcJ6Q0Q8BazO1pV1ErAttQPM/pJWSlojaaGk0oHLzMw6RxHRvS+TXgPOiogrCulPA9dFxLkly7kLGB0RBxfSzwT+ROrj2RWYDhwMHBYRDwxR1jRgGsDYsWMPnj17dnM7NYRVq1ax4447dqSskcztVNsjz6zYaHnsKHj+lQ3LB+45uss16g0+nsrpZDtNnDjxwYjoq7Wu56aKkbQH6XLa2cV1EfGNQt5bgV8B55JuCthERFwDXAPQ19cX/f39HannwMAAnSprJHM71Ta1xlQxlz2y4b/r8o/0d7lGvcHHUzndaqduXyIbBGr96TUmW1fGhwEBP2yUMSJWA7cCf1O2gmZm1hndDjBLKfS1SNoL2IFC30wdpwILI+J3JfNH9jIzsy7qdoCZAxwjaadc2hTgFWBBo40ljQMOpeTdY5JGAccDDzZbUTMza0+3A8xVwBrgRklHZh3sM4DL87cuS1om6bs1tj8VWAv8a3GFpNGS7pH0SUlHSJoCzAfeBlxcwb6YmVkdXe3kj4hBSUcAVwI3Ay8BXycFmWK9ak0fcypwZ0S8WGPdGuAPpBkBdgNeBe4HDo+IRR3ZATMzK63rd5FFxKPABxrkGTdE+l/V2eZVYHJblTMzs47xbMpmZlYJBxgzM6uEA4yZmVXCAcbMzCrhAGNmZpVwgDEzs0o4wJiZWSUcYMzMrBIOMGZmVommAoykWtO3mJmZbaLZM5hnJF0q6S8qqY2ZmY0YzQaYq4CTgV9K+qmkaZJ2rqBeZmbW45oKMBExIyL2AY4CHgMuB56T9H1JR1ZRQTMz600tdfJHxF0R8ffA7sDngHcDcyUtlzRD0ts6WUkzM+s97d5F1ge8n/QY5EHgHuDjwDJJH22zbDMz62FNBxhJe0u6QNITwJ3AHsA/AG+LiP8B7A1cDXytozU1M7Oe0tQDxyTNB94HPAPMBGZGxJP5PBHxuqQfAGd2rJZmZtZzmn2i5QvAJGBeRESdfIuBd7RcKzMz63nNXiL7FnBfreAiaUdJ7weIiNeKZzZmZrZlaTbAzAf2H2Ldu7P1ZmZmTQcY1Vm3I7C6jbqYmdkI0rAPJrvs1Z9L+rikYwvZtgeOBx7pXNXMzKyXlenkfy9pMCVAAKcAawt5/gQsBc7qXNXMzKyXNQwwEfE1sjEtkn4LfCgiFlddMTMz621N3aYcEb712MzMSinTBzMJWBgRK7PPdUXErR2pmZmZ9bQyZzA/Bg4FHsg+B0PfTRaAH0pmZmalAsw7gOdyn83MzBoq08n/ZK3PZmZm9ZTpg9mhmQIjwoMtzcys1CWyVaS+lbLcB2NmZqUCzD/QXIAxMzMr1QdzbRfqYWZmI0y7j0w2MzOrqUwn/wPA1Ih4VNLPaHC5LCIO6VTlzMysd5Xpg/kV8Erus/tjzMysoTJ9MB/LfZ5aaW3MzGzEaLkPRsmukuo9hMzMzLZQTQcYSZMk3Qe8CvweeFXSfZKO73jtzMysZzUVYCR9EriZNPjyTNLDx87Mlv89W29mZtbc82CAc4GrI+J/FtKvknQV8EXg6o7UzMzMelqzl8jeCvzbEOtuAN7SqABJ+0u6U9JqSc9KukhS3ellJI2TFDVes2vkPUHSI5JelfSopCml9szMzDqq2TOY+cDhwLwa6w4H7q63saQxwB3Ao8AJwL7AZaRAd16J7/88cG9u+cVC+YeRAt23gTOAScAsSYMRcXuJ8s3MrEPKDLTcP7f4TeA7kt4K3AS8AOwGfAg4Dvh4g+I+BYwCJkfESmCepJ2BGZIuzdLqeSwiflJn/ZeAuyPijGx5vqQDgPMBBxgzsy4qcwbzSzYeXCngk9mr+HTL26g/m/JxwNxCIJkNXEI6A7q5RH1qkrQdMJF05pI3G5gpaXRErGi1fDMza06ZADOxg983HrgrnxART0lana1rFGBmSnoL6cxpFvDFiFg/y8C+wLbA0sI2S0iX4N4F/Ky96puZWVllRvIv6OD3jQFeqpE+mK0byhrgW6TLXCuBfuBsUlA5IVc2NcofLKzfiKRpwDSAsWPHMjAwUK/+pa1atapjZY1kbqfaph+4dqPlsaM2TnOb1ebjqZxutVOznfxvkLQVsH0xvYonWkbEc8Bnc0kDkp4Hvi3pLyPi4TbKvga4BqCvry/6+/vbqusbFRwYoFNljWRup9qmnnPLRsvTD1zLZY9s+O+6/CP9Xa5Rb/DxVE632qnZgZaSdLakZcBrwB9rvOoZBEbXSB/DhjONsn6UvR+cK5sa5Y8prDczsy5odhzMGcA5wHdJnfv/C7gIeBxYTnapqY6lpL6WN0jaC9iBTftOGonC+xOkoDe+kG88sC6ro5mZdUmzAeYTwAXApdnyTRFxIXAAKUDs12D7OcAxknbKpU0hPQ6g2b6ek7P3BwEiYg1pnM4phXxTgPt9B5mZWXc12wfzDmBxRLwu6TXgzQARsU7St4HvkM5whnIV6SzoRkmXAPsAM4DL87cuZ5fgFkTE6dnyDGAn0iDLlcD7gbOAGyPiF7nyv0zqn7mCNE5nUvY6tsn9NDOzNjV7BvMfwI7Z56eAv86tG0MaRDmkiBgEjiCNlbkZuBD4OumsKG8bNh5Ps5Q0TmYmcCtwGvC17D1f/kLSmc2RwFzgg8BpHsVvZtZ9zZ7B3Au8h/Qj/wPSCPy3AH8CPgPc2aiAiHgU+ECDPOMKy7NJAyYbioibSGcvZmY2jJoNMDOAPbPPF5MukU0lnbnMAz7XqYqZmVlvayrARMRjwGPZ5zWkZ8GcWUG9zMysx7Uz0PLtwB7AsxHxTOeqZGZmI0Erj0z+tKTfAU8CPwWekvS0pOJDyMzMbAvW7Ej+84ErSeNZjgf6svc5wDez9WZmZk1fIvsMcHFEfKmQfls2N9hnSCP7zcxsC9fsJbJRDP3UygXUmPzSzMy2TM0GmJuAyUOsOwn4cXvVMTOzkaLMI5Mn5RbnAJdKGsemj0w+APjnzlfRzMx6UZk+mB+z6aOR9wSOqZH3e6QnTZqZ2RauTIB5R+W1MDOzEafMI5Of7EZFzMxsZGl6JL+kbUgd+ocBbwH+E7iHNHX+2nrbmpnZlqOpACNpN+B24CDSEyyfByaQxr88LOnoiPhDpytpZma9p9nblC8H3gocGhH7RMSEiNgHeG+WfnmnK2hmZr2p2QAzCTg7Ih7IJ0bEz4AvkKaNMTMzazrAbAf8cYh1fwTe1F51zMxspGg2wPwEOFvSn+UTs+Wzs/VmZmZN30U2HZgP/E7S7aRO/t1Igy4F9He0dmZm1rOaOoOJiMXAfsA1wK7AUaQAcxWwX0Q83PEamplZTyp9BiNpW+AQ4LcRcU51VTIzs5GgmTOY14G7gPEV1cXMzEaQ0gEmItYBvwZ2r646ZmY2UjR7F9kXgfMlHVhFZczMbORo9i6y80gj9hdLeoZ0F1nkM0TEIR2qm5mZ9bBmA8wvs5eZmVldpQKMpFGkaWJ+CfweuCMinq+yYmZm1tvKPDJ5H+AOYFwueaWkD0fE7VVVzMzMeluZTv5LgXXA+4AdgAOAh4CrK6yXmZn1uDIBZgJwXkTcGxGvRsQS4JPAn0vao9rqmZlZryoTYPYAflNIe4I095jHxJiZWU1lx8FE4yxmZmYblL1Nea6ktTXS7yymR8Ru7VfLzMx6XZkAc2HltTAzsxGnYYCJCAcYMzNrWrNzkZmZmZXiAGNmZpVwgDEzs0o4wJiZWSUcYMzMrBIOMGZmVomuBxhJ+0u6U9JqSc9KukjS1g22eY+kmZKWZds9JukCSdsX8s2QFDVex1a7V2ZmVtTsA8faImkMaer/R4ETgH2By0iB7rw6m07J8l4C/Bo4CPhy9n5SIe8KoBhQlrRbdzMza05XAwzwKWAUMDkiVgLzJO0MzJB0aZZWy1cj4sXc8oCkV4GrJe0dEU/m1q2NiJ9UU30zMyur25fIjgPmFgLJbFLQOXyojQrBZb2Hsve3da56ZmbWKd0OMOOBpfmEiHgKWJ2ta8YE0oPQniikv1nSi5Jek/SQpMkt19bMzFqmiO7NxC/pNeCsiLiikP40cF1EnFuynN2BXwC3RsTUXPpHgd1IZzc7kR6MNgk4KSJuHKKsacA0gLFjxx48e/bsZnerplWrVrHjjjt2pKyRzO1U2yPPrNhoeewoeP6VDcsH7jm6yzXqDT6eyulkO02cOPHBiOirta7nAoykN5FuFHg7cHBEDNbJK+A+YFRE/FWjsvv6+mLRokWNspUyMDBAf39/R8oaydxOtY0755aNlqcfuJbLHtnQZbr8q8d3u0o9wcdTOZ1sJ0lDBphuXyIbBGr96TUmW1dXFjCuAw4AJtULLgCRoueNwEGNboU2M7PO6vZdZEsp9LVI2gvYgULfzBCuIN3efFRElMkP6WmcfiKnmVmXdfsMZg5wjKSdcmlTgFeABfU2lPQF4LPARyNiYZkvy854TgIejojXW6uymZm1ottnMFcBZwA3SroE2AeYAVyev3VZ0jJgQUScni2fBlwMXAs8I+nQXJlPRMQfsnwLgBtIZ0N/BnwCeC9wYrW7ZWZmRV0NMBExKOkI4ErgZuAl4OukIFOsV77P5OjsfWr2yvsYKfAALAP+EdiDdAvzz4HjI2JOJ+pvZmbldfsMhoh4FPhAgzzjCstT2TSw1Nru9DaqZmZmHeTZlM3MrBIOMGZmVgkHGDMzq4QDjJmZVcIBxszMKuEAY2ZmlXCAMTOzSjjAmJlZJRxgzMysEg4wZmZWCQcYMzOrhAOMmZlVwgHGzMwq4QBjZmaVcIAxM7NKOMCYmVklHGDMzKwSDjBmZlYJBxgzM6uEA4yZmVXCAcbMzCrhAGNmZpVwgDEzs0o4wJiZWSUcYMzMrBLbDHcFzGz4jTvnloZ5ln/1+C7UxEYSn8GYmVklHGDMzKwSDjBmZlYJBxgzM6uEA4yZmVXCAcbMzCrhAGNmZpVwgDEzs0o4wJiZWSUcYMzMrBIOMGZmVgkHGDMzq4QDjJmZVcKzKZs10GimYc8ybFabz2DMzKwSDjBmZlaJrgcYSftLulPSaknPSrpI0tYlthstaaakQUkrJH1f0ltr5DtB0iOSXpX0qKQp1eyJmZnV09U+GEljgDuAR4ETgH2By0iB7rwGm18PvAv4OLAOuAS4CXhfrvzDgBuAbwNnAJOAWZIGI+L2ju6MmbXMT9DcMnS7k/9TwChgckSsBOZJ2hmYIenSLG0TkiYARwOHR8TdWdozwE8lHRkRd2RZvwTcHRFnZMvzJR0AnA84wIwQ/nEy6w3dDjDHAXMLgWQ26WzkcODmOts9vz64AETEA5J+m627Q9J2wETSmUvebGCmpNERsaJD+2FmPcB3AA6vbgeY8cBd+YSIeErS6mzdUAFmPLC0RvqSbB2ky23b1si3hHQJ7l3Az1qrdnPGnXML0w9cy9QSf2kPpdaB36n/LFWdAfjMwrYUnTzWW/l/Xeb765l+4Fr62yqhHEVEF74m+zLpNeCsiLiikP40cF1EnDvEdvOAlyPixEL694B9IuK/SvpbYCHw1xGxOJfnncCvgWNq9cNImgZMyxbfDTzW8g5ubBfgxQ6VNZK5ncpxO5Xjdiqnk+20d0TsWmvFFj/QMiKuAa7pdLmSFkVEX6fLHWncTuW4ncpxO5XTrXbq9m3Kg8DoGuljsnXtbLf+vZhvTGG9mZl1QbcDzFI29JkAIGkvYAdq97EMuV0m3zfzBPBajXzjSbc1P95Cfc3MrEXdDjBzgGMk7ZRLmwK8AixosN3u2TgXACT1Aftk64iINcB84JTCtlOA+4fhDrKOX3YbodxO5bidynE7ldOVdup2J/8Y0iDLX5JuTd4HuBy4IiLOy+VbBiyIiNNzaXOB/YDPs2Gg5QsRURxoOQBcSRqEOSnLf6wHWpqZdVdXz2AiYhA4AtiadEvyhcDXgQsKWbfJ8uRNIZ3l/AtwHfAg8KFC+QuBk4EjgbnAB4HTHFzMzLqvq2cwZma25fBsyg14cs7GWmkjSe/J2mdZtt1jki6QtH0h3wxJUeN1bLV71XktttO4IfZ/do28PX8sQcvtNNRxEpK+kMt37RB5at1EtFmT9E5JV0v6haTXJQ2U3K5rv01b/DiYejw5Z2NttNGULO8lpIGwBwFfzt5PKuRdARQDypJ2695NbR5LkPoS780tbzRIbiQcS9BWO30HuK2QdiJwNtmNQDlLgY8V0pa3VuNhdQDp3/knpFlMyureb1NE+DXEC/gCafzMzrm0fwZW59NqbDcBCOD9ubRDsrQjc2lzgbsK294KLBzufe9CG+1SI21a1kZ759JmAC8O934OYzuNy9rk7xqU3/PHUjvtNERZtwBLCmnXAouGez871FZb5T7/CBgosU1Xf5t8iay+oSbnHEWanLPedptMzgmsn5yT3OSc1xe2nQ1MkFRrYOnmqKU2ioha01Q8lL2/rXPV22y0eiw1NIKOJehQO2WXfI4CZnW2epuPiFjXwmZd/W1ygKlvk0k2I+Ip0l9T9a7Zdmpyzl7QahvVMoF0yv5EIf3Nkl6U9JqkhyRNbrm2w6fddpqZXWd/TtLlkkbl1o2UYwk6dzydRGqTWgFmf0krJa2RtFBSWwG+x3T1t8kBpr4xwEs10gfZMAVNq9utfy/mGyys39y12kYbkbQ76Rr7/42IF3KrlpEukZxC+tF4FrihB4NMq+20BvgWcDrpFv+rgU+T/prMl02N8nvtWIIOHU/AqcDPI+LXhfSHgOnAfwM+QhoOMU/SIS3UtRd19bfJnfw27CS9iXQ6vgr4p/y6iPheIe/NwH2kh8jd2K06DpeIeA74bC5pQNLzwLcl/WVEPDxMVdtsSdqDdDnt7OK6iPhGIe+twK+Ac0k3BVgH+QymPk/O2VirbQSAJJEGzh4ATIo0GHdIkXobbwQOKnO7+GakrXYq+FH2fnCubGqU32vHEnSmnT4MCPhho4wRsZrUef03ZSvY47r62+QAU58n52ys1TZa7wrS7agnRESZ/JDueOm1EcLttlNeFN5HyrEEnWmnU0l3O/2uZP5ePJ5a1dXfJgeY+rakyTlb1WobkQ2A+yzw0UjT/DSUnfGcBDwcEa+3VuVh0XI71XBy9v4gjKhjCdpsJ0njgEMpefdYdrPE8WRtuQXo7m/TcN/LvTm/SKeEzwHzSPObTSP1E3ylkG8Z8N1C2lzgN8Bk0rXdx4B7CnkOA9aS/orvBy4l/YVw9HDve9VtBJxG+qtxJukHIf/aNZdvAWmg19Gkueduzdrog8O9711qpxmkgYaTs+0uIv3Y3jDSjqV22imXfg7pr+9a46xGA/cAnyTdMDGFNEhxDdA33PveQlvtQPpj42TgflJf0vrlHYZqp27+Ng17I23uL2B/4K7sP/VzpNHmWxfyLAeuLaS9OfvxfAlYCfxgiIP+RNLs0mtIp6inDvc+d6ONSAPeYojX1Fy+72b/GV4BXs5+II4b7n3uYjudCiwizWbwp+wH4yJgu5F4LLXaTrn0xcBtQ5S7Pan/7ndZG60gjf4/dLj3ucV2Glfn/9C4odqpm79NnuzSzMwq4T4YMzOrhAOMmZlVwgHGzMwq4QBjZmaVcIAxM7NKOMCYmVklHGDMzKwSDjBmZlaJ/w97k4So4+RPgQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZgAAAEPCAYAAAB/WNKuAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3deZxcVZ338c+XZSRsISwJoEiMI+YBdZS0LDOMBAGVZQyCGAYdHwSSOIODMy9ABFGCqM8Asjg4PiSoIMqAC4qjLCFAOogLkAR4kCUQNKyyB0JIgAR+zx/nNlxuqqtvVdetrup836/XfVXXuefeOqe6u351z3YVEZiZmbXaWkNdADMzG54cYMzMrBIOMGZmVgkHGDMzq4QDjJmZVcIBxszMKuEAYy0labqk6Gf7VMlzbJedZ5NC+mHZeTaspvTlyjHIc/5MUm+JfOtI+jdJt0taIWmJpCsl7dbk63bKe3pY4W/iCUmzJO1Y4tiJ2THvakdZbfAcYKwKzwG71tiuLnn8dsDJQPGD/YrsPMtbU8ymy1EpSWsDlwPfAP4H2Bc4DHgF6JV0aBOn7ZT3tM8Hs9edBmwBzJG09QDHLMiOub/islmLrDPUBbBhaVVE/KHVJ42IJ4EnW33eDvSvwH7APhGRD8q/lHQpMFPS3Ih4ZLAvNITv6S0RsQxA0jzgAeCTwBnFjJIEvCkilgIt/7uy6vgKxoaEpBMkLZL0oqTHJV0taUtJE4FfZdn+nDWJLM6OeUNzjqSx2fNDJF0gaamkh/ua4iR9QdKjkp6UdJqktXKvP17SpZIekrRc0p1Zk9Ra2f5+y5Htf2t2/DPZ8bMkvbNQx22yZq0VkhZLOrLk2/N5YE4huPT5ErAecETudRZL+qakL0t6TNIySRdLGjlQXWo1kUnaXNIPJD2d1a1XUk+hbn2v+e/Ze74kez8avtqLiIdIQW5sdu7pkp6StJukW4AXgYNrNZFJWjv7W7pX0ktZWS4slHWSpHnZ39pjkk6XtG6j5bTG+QrGKiFptb+tiFiV7fs0cCJwPHAnsBmpyWQDUjPIscA3gQOBvwAvDfBypwEXAwcBhwM/kPQ+YNvs+QTga8CtwKXZMW8GFmbHPQ+8FzgFGAH8n3rlkLQpcCPwNPBZUvPSF4FrJW0XESuyb92/BDYnBYMXs/NvCtxX533bhvRBe3at/RFxv6Q7gA8Udv0jsAiYAmwFnA58Fzi4Xl36cTnw19kxTwHHkZqw3hcRi3L5PgH8P2Aq8BbgLFKz3r/UOfdqJG1Eel8eyyWvD/wgq8e9wKNZvYpmAJ/O8s3NznNQ7tyfAC7J8p0IvJ30+10rq59VKSK8eWvZBkwHop9tbJbn28Bldc6xfz5/Lv2wLH3D7PnY7PkFuTwbAytJH+Jr59JvBn7cz+uJ9GXrROBPJcpxKim4bJpLG0Xqezoqe75vduzOuTzbAquA3jp13yU7blKdPJcDd+eeLwae6XtfsrRPAq8C/6vB9/Qj2fPdc3k2IF1hzCi85v3AOrm0c4DHBvj76Hu9kdl7vg3w4+x9eW/hb2hS4diJWfq7sufjs+dH1/m9PpD/+8jSDwdWAJsN9f/LcN98BWNVeA7Yq0b6o9njbcARkk4hdTLPj4hXBvF61/X9EBFLJT0JzC2ccxHw1r4nktYDTiB9EL8VWDe3b53Irrb6sRcwG1iau1J7HpgP9DUl7QQ8HhE35cr2gKT5TdSvjNmR9WlkfkH6gH0/cHcD59kJeCIi5vYlRMQLkn4NFEewzSm8T3cBoyWtGxErB3idZ3M/PwUcHhG35dICuGqAc+yRPV7Yz/7tSL/bnxSuqK8nNTO+i3TVYxVxgLEqrIqIeXX2fx/YiNS08hXgaUnnASc3GWieLTx/uZ+09XLPTwOOJDVbLcjyTwJOyvIto3+bk640JtfY1xfstgSeqLH/CVLd+9PXcb9tnTzb5vLlz/uaiFguaRm1m5Xq2ap4rszjpOanvFrvsYA3ka4i6/kAqWnxKeChiHi1sH9JRLw8wDk2A16I1Plfy+bZ45X97N9mgPPbIDnAWNtlHyZnA2dnfQ6fBL4OPAyc16ZiHAycGxGn9yVI2q/ksc+Qhg+fWmPf89njY8DoGvtHk5pnaoqIh7IO+I8C/1ncL+ltpG/exdceXci3PrAhqb+lEX8pniszhlTvVrm1cMVVVOY+Ik8DG0jauJ8g01feqaT+t6I/l3gNGwSPIrMhFREPRcR/kJqwts+S+765rlf7qJYYQa6jW2nuySGFPP2V4zpgB+DOiJhX2BZmeW4BxkjaOfcabwUGnFAIfAvYU9KHauz7elbu7xXS99YbJ0t+jPQh3XclWfY9vYnUzPXaIIIsWO1HGtjQSa7PHj/dz/6FpCu9sTV+T/Mi4un2FHPN5SsYq8I6knapkf5QRDwiaQbp2+UfSP01ewDvII0qg/TBADBNad7H8oi4o8VlnA0cJWlRVpajSE07ef2V4yzgU8D1ks4lfYiNAXYHboyIS0jNMrcDP5V0PCkonELt5qeic0n9PL+Q9E2gl9SsdgSps/6fYvU5MCuAKySdQWrmOgP4RUTcNUBd3iAiZkn6HfBjSV8kXSUcSwrIq81RGUoRsVDSTOBMSaOBG0gTST8eEYdExKuSjgF+KGljUp/Oy8A44IAsX7snmK5ZhnqUgbfhtVF/FNlJWZ7DgN+SPtiXk4a6HlE4zzGkEUCrgMW542qNItu/cOxi4JuFtAuBebnnY0gd4UtJ/Qunk4b4vnb+/sqRpW8NXJAd+1L2mj8CdsjleStp9YIV2TmmAT+jziiy3LHrAP+evTcrgCWkD8jdauRdDJyZvfePAy+QhuZu0uh7mqVtAVyUveYKUkf4+0u8x6udq0ZZy+SZDjxVI30iuVFkWdraZKP/SMHjYeD7heP2AX6TvS9LSYNMvkZuBJy3ajZlv4C2kfTXpHH1u5KaGX4TERNLHDeSNAzyAFLT3q9JwxOfLuSbRPrjeQfpj+6UiPhxK+tg1kmyPpufRYTndVhHGYo+mB1IcwQWkiZQlfUT0jeYI0nfgt5Pmg/wGqWFAC8D5pC+tVwBXNJPW7aZmVVoKK5g1opsSKKknwGbD3QFI2lX4HekyV83ZGk7kTok946Ia7O0WcC6EfHB3LFXAhtHRFOr0Jp1Ol/BWKdq+xVMrD7evYx9SJPWbsid52bSMMN9ACS9idRZ/JPCsZcCu/aty2Q23ETEWAcX60TdMkx5PHBPjfS7s32Q1hhat0a+u0n13K6y0pmZ2Wq6ZZjyKFafNQxplMu4XB5q5FtS2P8GkqaSJmIxYsSICdtsU25y76uvvspaa3VLfK5tONQBXI+N7k1dmc9vN/Tfodb030WnaUc97r333qciYota+7olwFQmImYCMwF6enpi3rx6K5y8rre3l4kTJ1ZYsuoNhzqA64GUHhcurJ+vDdb430WHaUc9JD3Q375uCdFLSKuvFo3i9SuUvsdivlGF/WZm1gbdEmDu4fW+lrx838z9pAX2ivnGk5Ytb2RItJmZDVK3BJirgC2zeS4AZHfYG5ftIyJeIs1/Obhw7GTg9xHxXJvKamZmDEEfTLZw3r7Z0zcDG0v6ePb8ykjLjC8i3c/jCICI+L2ka4CLJB1LuiI5jbTu07W5058K9Eo6hzQJc99s+0jlFTMzszcYik7+0cBPC2l9z99GWuNoHdIaQ3mTSUu8f5/cUjH5DBFxYxasvgb8M2mezKERcU0Ly2/WWdo8WdqsrLYHmIhYTLopUb08Y2ukPQt8JtvqHXs5hSVkzMys/bqlD8bMzLqMA4xZt5swIW1mHWaNn2hp1vUWLBjqEpjV5CsYMzOrhAOMmZlVwgHGzMwq4QBjZmaVcIAxM7NKeBSZWbebMmWoS2BWkwOMWbebOXOoS2BWk5vIzMysEg4wZt1u/vy0mXUYN5GZdbuenvToVZWtw/gKxszMKuEAY2ZmlXCAMTOzSjjAmJlZJRxgzMysEg4wZmZWCQ9TNut28+YNdQnManKAMet2vl2ydSg3kZmZWSUcYMy63dSpaTPrMA4wZt3u/PPTZtZhHGDMzKwSDjBmZlYJBxgzM6uEA4yZmVXCAcbMzCrhiZZm3W7HHYe6BGY1OcCYdTvfLtk6lJvIzMysEg4wZmZWCQcYs24npc2swzjAmJlZJRxgzMysEg4wZmZWCQ9TNusQdzzyHId98YqGj1vc+qKYtYSvYMzMrBIOMGZmVgk3kZl1uxkzhroEZjU5wJh1O98u2TpU25vIJG0v6TpJyyU9KumrktYe4JjpkqKf7YRcvgv7yTO++pqZmVleW69gJI0CrgXuAiYBbwfOJAW6k+oc+l3g6kLaAcDxwFWF9HuAzxTSFjdXYrMuMHNmevSVjHWYdjeRfRYYARwYEUuB2ZI2BqZLOj1LW01EPAw8nE+T9GXgnoi4rZD9hYj4QwVlN+tM06alRwcY6zDtbiLbB5hVCCSXkoLO7mVPImkzYG/gktYWz8zMWqXdAWY8qQnrNRHxILA821fWQcC61A4w20taKuklSTdKKh24zMysddrdRDYKeLZG+pJsX1mHAAsi4r5C+q3ATaQ+ni2AY0jNcLtFxM21TiRpKjAVYMyYMfT29pYqwLJly0rn7VTDoQ4wfOoxZgQc8+5VTR/fCe/BcPlduB6t0XXDlCVtRWpOO764LyK+Vch7JXAncCJpUMBqImImMBOgp6cnJk6cWKocvb29lM3bqYZDHWD41OPci3/JmXc0/i/5r9ljJ7wHw+V34Xq0RrubyJYAI2ukj8r2lfEJQMCPB8oYEcuBKwHftNzMrM3aHWDuodDXImkbYH0KfTN1HALcGBEPlcwf2WZmZm3U7gBzFfBhSRvl0iYDK4C5Ax0saSywCyVHj0kaAewHzG+0oGZdIyJtZh2m3QHmPOAl4OeS9so62KcDZ+WHLktaJOl7NY4/BFgF/LS4Q9JISb+RNE3SnpImA3OArYFvVFAXMzOro62d/BGxRNKewLeBX5FGlJ1NCjLFctVaPuYQ4LqIeKrGvpeAJ0krAowGXgR+D+weEfNaUgEzMyut7aPIIuIu4IMD5BnbT/p76xzzInDgoApn1o0mTEiP890SbJ2l64Ypm1nBggVDXQKzmnzDMTMzq4QDjJmZVcIBxszMKuEAY2ZmlXCAMTOzSngUmVm3mzJlqEtgVpMDjFm367tlslmHcROZmZlVoqEAI6nW8i1mNpTmz/csfutIjTaRPSLpIuCCiLi7igKZWYN6etKjV1S2DtNoE9l5wMeBP0q6SdJUSRtXUC4zM+tyDQWYiJgeEeOAvYGFwFnAXyRdLGmvKgpoZmbdqalO/oi4PiI+DWxJuiX4O4FZkhZLmi5p61YW0szMus9gR5H1AB8g3QZ5CfAb4EhgkaRPDfLcZmbWxRoOMJK2lXSypPuB64CtgMOBrSPin4BtgRnAGS0tqZmZdZWGRpFJmgP8PfAIcAFpNNkD+TwR8Yqk/wY+37JSmplZ12l0mPITwL7A7Ii6YyJvA97WdKnMrLx5viO4daZGA8x/AQtqBRdJGwI7RsQNEbESeGC1o82s9fpumWzWYRrtg5kDbN/Pvndm+83MzBoOMKqzb0Ng+SDKYmbNmDo1bWYdZsAmMkkfACbmko6U9JFCtvWA/YA7Wlc0Myvl/PPTo1dVtg5Tpg9mZ9JkSoAADgZWFfK8DNwDHNe6opmZWTcbMMBExBlkc1ok/Rn4WETcVnXBzMysuzU0iiwiPPTYzMxKKdMHsy9wY0QszX6uKyKubEnJzMysq5W5gvk1sAtwc/Zz0P9osgB8UzIzMysVYN4G/CX3s5l1kh13HOoSmNVUppP/gVo/m1mH8O2SrUOV6YNZv5ETRoQnW5qZWakmsmWkvpWy3AdjZmalAszhNBZgzKydlI25qbvAuVn7lemDubAN5TAzs2FmsLdMNjMzq6lMJ//NwGERcZekWxiguSwidmpV4czMrHuV6YO5E1iR+9kNvWZmNqAyfTCfyf18WKWlMTOzYaPpPhglW0iqdxMyMzNbQzW0mjK8tvjlScCE7PhVkuYDX4+IK1pcPjMbyIwZQ10Cs5oaCjCSpgHfAa4DPg88AYwGDgT+R9K/RIT/2s3aybdLtg7V6BXMicCMiPiXQvp5ks4DvgQ4wJiZWcN9MJsBv+hn32XApgOdQNL2kq6TtFzSo5K+Kqnu8jKSxkqKGtulNfJOknSHpBcl3SVpcqmamXWrmTPTZtZhGr2CmQPsDsyusW934IZ6B0saBVwL3AVMAt4OnEkKdCeVeP1jgd/mnj9VOP9upED3HeBoYF/gEklLIuKaEuc36z7TpqVHN5VZhykz0XL73NP/BL4raTPgcl7vg/kYsA9w5ACn+ywwAjgwIpYCsyVtDEyXdHqWVs/CiPhDnf1fBm6IiKOz53Mk7QB8BXCAMTNrozJXMH/kjZMrBUzLtuLdLa+m/mrK+wCzCoHkUuA00hXQr0qUpyZJbwL2IF255F0KXCBpZEQ81+z5zcysMWUCzB4tfL3xwPX5hIh4UNLybN9AAeYCSZuSrpwuAb4UEX2rDLwdWBe4p3DM3aQmuO2AWwZXfDMzK6vMTP65LXy9UcCzNdKXZPv68xLwX6RmrqXAROB4UlCZlDs3Nc6/pLD/DSRNBaYCjBkzht7e3nrlf82yZctK5+1Uw6EOMHzqMWYEHPPuVU0f3wnvwXD5XbgerdHwRMs+ktYC1iumV3FHy4j4C/C5XFKvpMeB70j6m4i4fRDnngnMBOjp6YmJEyeWOq63t5eyeTvVcKgDDJ96nHvxLznzjsb/Jf81e+yE92C4/C5cj9ZoaJhytjzM8ZIWASuB52ts9SwBRtZIH8XrVxpl/Sx7nJA7NzXOP6qw38zM2qDReTBHA18Evkfq3P868FXgXmAxWVNTHfeQ+lpeI2kbYH1W7zsZSBQe7ycFvfGFfOOBV7Mymg0/Eb6bpXWkRgPMFOBk4PTs+eURcQqwAylAvGOA468CPixpo1zaZNLtABrt6/l49jgfICJeIs3TObiQbzLwe48gMzNrr0YbfN8G3BYRr0haCWwCEBGvSvoO8F3SFU5/ziNdBf1c0mnAOGA6cFZ+6HLWBDc3Io7Ink8HNiJNslwKfAA4Dvh5RPy/3PlPJfXPnEOap7Nvtn2kwXqamdkgNXoF8zSwYfbzg8D7cvtGkSZR9isilgB7kubK/Ao4BTibdFWUtw5vnE9zD2mezAXAlcChwBnZY/78N5KubPYCZgEfBQ71LH4b1iZMSJtZh2n0Cua3wPtJH/L/TZqBvynwMnAUaZXluiLiLuCDA+QZW3h+KWnC5IAi4nLS1YvZmmHBgqEugVlNjQaY6cCbs5+/QWoiO4x05TKb10dMmpnZGq6hABMRC4GF2c8vke4J8/kKymVmZl1uMBMt3wJsBTwaEY+0rkhmZjYcNNrJj6R/lvQQ8ABwE/CgpIclFW9CZmZma7BGZ/J/Bfg2aT7LfkBP9ngV8J/ZfjMzs4abyI4CvhERXy6kX52tDXYUaWa/mbXLlClDXQKzmhoNMCPo/66Vc/EoMrP28+2SrUM12gdzOXBgP/sOAn49uOKYmdlwUeaWyfvmnl4FnC5pLKvfMnkH4AutL6KZ1TV/fnr0bH7rMGWayH7N6rdGfjPw4Rp5f0S606SZtUtPT3r0isrWYcoEmLdVXgozMxt2ytwy+YF2FMTMzIaXhmfyS1qH1KG/G7Ap8AzwG9LS+c3fUNzMzIaVhgKMpNHANcB7SHewfBzYlTT/5XZJH4qIJ1tdSDMz6z6NDlM+C9gM2CUixkXErhExDtg5Sz+r1QU0M7Pu1GiA2Rc4PiJuzidGxC3ACaRlY8zMzBrug3kT8Hw/+54H/mpwxTGzhs2bN9QlMKup0QDzB+B4SddHxAt9iZI2AI7P9ptZO3mCpXWoRgPMMcAc4CFJ15A6+UeTJl0KmNjS0pmZWddqqA8mIm4D3gHMBLYA9iYFmPOAd0TE7S0voZnVN3Vq2sw6TOkrGEnrAjsBf46IL1ZXJDNryPnnp0evqmwdppErmFeA64HxFZXFzMyGkdIBJiJeBe4DtqyuOGZmNlw0Og/mS8BXJL27isKYmdnw0egospNIM/Zvk/QIaRTZG9YIj4idWlQ2MzPrYo0GmD9mm5mZWV2lAoykEaRlYv4IPAZcGxGPV1kwMytpxx2HugRmNZW5ZfI44FpgbC55qaRPRMQ1VRXMzErqu2WyWYcp08l/OvAq8PfA+sAOwK3AjArLZWZmXa5MgNkVOCkifhsRL0bE3cA04K2Stqq2eGZm1q3KBJitgD8V0u4nrT3mOTFmQ01Km1mHKTsPJgbOYmZm9rqyw5RnSVpVI/26YnpEjB58sczMrNuVCTCnVF4KMzMbdgYMMBHhAGNmZg1rdC0yMzOzUhxgzMysEo2uRWZmnWaG5zxbZ3KAMet2vl2ydSg3kZmZWSUcYMy63cyZaTPrMG0PMJK2l3SdpOWSHpX0VUlrD3DM+yVdIGlRdtxCSSdLWq+Qb7qkqLF9pNpamQ2hadPSZtZh2toHI2kUaen/u4BJwNuBM0mB7qQ6h07O8p4G3Ae8Bzg1ezyokPc5oBhQ7h5s2c3MrDHt7uT/LDACODAilgKzJW0MTJd0epZWy39ExFO5572SXgRmSNo2Ih7I7VsVEX+opvhmZlZWu5vI9gFmFQLJpaSgs3t/BxWCS59bs8etW1c8MzNrlXYHmPHAPfmEiHgQWJ7ta8SupBuh3V9I30TSU5JWSrpV0oFNl9bMzJrW7iayUcCzNdKXZPtKkbQlqc/mhxHxRG7XIuALpKubjUg3RrtM0kER8fN+zjUVmAowZswYent7S5Vh2bJlpfN2quFQBxg+9RgzAo55d61Fy8vphPdguPwuXI/WUET7bvUiaSVwXEScU0h/GLgoIk4scY6/Ig0UeAswISKW1Mkr4HfAiIh470Dn7unpiXnz5g2UDUj/zBMnTiyVt1MNhzrA8KnHuRf/kjPvaPw73+LT9k8/tPF/uT/D5XfhepQnaX5E9NTa1+4msiXAyBrpo7J9dWUB4yJgB2DfesEFIFL0/DnwnoGGQpt1rYiOCC5mRe1uIruHQl+LpG2A9Sn0zfTjHNLw5r0jokx+SHfj9H+fmVmbtfsK5irgw5I2yqVNBlYAc+sdKOkE4HPApyLixjIvll3xHATcHhGvNFdkMzNrRruvYM4DjgZ+Luk0YBwwHTgrP3RZ0iJgbkQckT0/FPgGcCHwiKRdcue8PyKezPLNBS4jXQ1tAEwBdgYOqLZaZkNowoT0OH/+0JbDrKCtASYilkjaE/g28CvSiLKzSUGmWK58n8mHssfDsi3vM6TAA2kU2b8BW5GGMC8A9ouIq1pRfrOOtGDBUJfArKa2L9cfEXcBHxwgz9jC88NYPbDUOu6IQRTNzMxayKspm5lZJRxgzMysEg4wZmZWCQcYMzOrRNs7+c2sxaZMGeoSmNXkAGPW7Xy7ZOtQbiIzM7NKOMCYdbv58z2L3zqSm8jMul1PtlK6V1S2DuMrGDMzq4QDjJmZVcIBxszMKuEAY2ZmlXCAMTOzSjjAmJlZJTxM2azbzZs31CUwq8kBxqzb9d0y2azDuInMzMwq4QBj1u2mTk2bWYdxgDHrduefnzazDuMAY2ZmlXCAMTOzSjjAmJlZJRxgzMysEg4wZmZWCU+0NOt2O+7Y1GFjv3hF0y+5+D/2a/pYW3M4wJh1O98u2TqUm8jMzKwSDjBmZlYJBxizbielzazDOMCYmVklHGDMzKwSDjBmZlYJBxgzM6uEA4yZmVXCAcbMzCrhmfxm3W7GjKEugVlNDjBm3c63S7YO5SYyMzOrhK9gzGroqpWGZ85Mj76SsQ7jAGPW7aZNS48OMNZh3ERmZmaVaHuAkbS9pOskLZf0qKSvSlq7xHEjJV0gaYmk5yRdLGmzGvkmSbpD0ouS7pI0uZqamJlZPW1tIpM0CrgWuAuYBLwdOJMU6E4a4PCfANsBRwKvAqcBlwN/nzv/bsBlwHeAo4F9gUskLYmIa1paGTMrrav6tKxl2t0H81lgBHBgRCwFZkvaGJgu6fQsbTWSdgU+BOweETdkaY8AN0naKyKuzbJ+GbghIo7Ons+RtAPwFcABpsv19yF1zLtXcdgAH2D+kDJrv3YHmH2AWYVAcinpamR34Fd1jnu8L7gARMTNkv6c7btW0puAPUhXLnmXAhdIGhkRz7WoHmbWYXyV1HnaHWDGA9fnEyLiQUnLs339BZjxwD010u/O9kFqblu3Rr67SU1w2wG3NFfsxgzmD70K+X+efNnKfPPv7zxF/ue24apVf9tD8blQ73+8Hf93iojKX+S1F5NWAsdFxDmF9IeBiyLixH6Omw28EBEHFNJ/BIyLiL+V9HfAjcD7IuK2XJ6/Bu4DPlyrH0bSVKBvfOc7gYUlq7M58FTJvJ1qONQBXI9OMhzqAK5HI7aNiC1q7Vjj58FExExgZqPHSZoXET0VFKlthkMdwPXoJMOhDuB6tEq7hykvAUbWSB+V7RvMcX2PxXyjCvvNzKwN2h1g7uH1PhMAJG0DrE/tPpZ+j8vk+2buB1bWyDeeNKz53ibKa2ZmTWp3gLkK+LCkjXJpk4EVwNwBjtsym+cCgKQeYFy2j4h4CZgDHFw4djLw+wpGkDXcrNaBhkMdwPXoJMOhDuB6tES7O/lHkSZZ/pE0NHkccBZwTkSclMu3CJgbEUfk0mYB7wCO5fWJlk9ERHGiZS/wbdIkzH2z/B/xREszs/Zq6xVMRCwB9gTWJg1JPgU4Gzi5kHWdLE/eZNJVzveBi4D5wMcK578R+DiwFzAL+ChwqIOLmVn7tfUKxszM1hxeTbkkSVMk3Zctojlf0p4NHv8+Sa9IGtKx9c3UQ9I0SbMlPZ4tNPpbSR9qQ1krXRi1XZqph6T3Z3VYlB23UNLJktZrV7lrlKmp30fu+LUkzZMUkvavsqx1ytB0HSQdKOkWSSskPS3pakkbVF3mfsrS7P9Gj6RrJD2TbddK2rmygkaEtwE24B+BV0hrne1BaqJbAbyr5PECfgs8BjzVbfUAHiR1Fh4A7A38gNQP9tEKyzoKeJS0OOrepHXsXgC+VuLYWcCfgYNIzaj3Ar8Zove8qXoA33XcRS4AAAWASURBVARuAKYAE0lLID0HXNZN9SicY2r2PxDA/t1UB9Iiuy8CX81+Hx8DzgVGdks9gG2AZ0mrqeyXbb3AUtJkydaXtd1vTjdupNn93889Xwu4A/hRyeP/CVgEfGOIA0xT9QA2r5H2O2BOhWU9gTR3aeNc2heA5fm0Gsftmn2AfSCXtlOWttcQvOfN1qPWez41q8e23VKPXN5RwJPAEUMYYJr+XQDPA1PaXeYW1+OzpC+YI3Npo7K0f66irG4iG4CkcaR1zH7SlxYRrwI/JS20OdDxG5FGvB0LvFxRMQc0mHpERK1mvVuBrVtZxoL+FkYdQVoYtd5xqy2MSrqiGfD3VYGm6lHnPYdq3/f+NPv76HMq6Sr+ugrKVlazdfhE9viDqgrWoGbrsS6winS102dZlqZWFxLcB1NG38TNWotobiqp5ho8OV8B7o6Iy1tessYMth5Fu1Lt5NXVFjiNiAdJ39JqTbrt97hMfmHUdmq2HrXsSmqavL81RWtI0/WQ9B7gcNKXrKHUbB12Jl39HyHpYUkrJd0k6W+rK2pdzdbjsizPmZJGSxpNGsW7hPRFs+UcYAbWt9TMs4X0JYX9q5H0TuAo4N8qKFejmq5HkaTDgfeR5jBVZRSrlxVSeeuVtdnjqtKS8kjaknRTvh9GxBMtKlsjBlOPc4FvR8SilpeqMc3WYUvSQrgnAccD/0C6Crha0phWF7KEpuoREY+S+l4PAh7PtgNJCwE/WUE518zFLiWNBLYaKF9E1Fu+poxvARdGxB2DPE9NbaxH/jUnkD4wvhURc1p1XuufpL8iNW0uA/59iIvTEEmHkD6c/2GoyzIIAjYEDo6IqwEk/Q54APgcadBMx5O0FelKZT5p0AKkL8BXSPrb7CqopdbIAENaTub8EvnEGxfRzH9rqLuIpqR9gL8DPidpkyx5vbRLmwArIi1vMxiV1+MNJ0n9OFeQ2tGPKV/MpgxmYdRazX0DHVeVZusBpD8W0mi/HYC/izRZeSg0XA9J6wJnkPog18r+7jfOdm8gaaOIeL6KwvZjMH9TQRpxBUBELJU0H9i+lQUsqdl6HEfqh/l4RKwEkHQ96XYmx7L6zRoHbY1sIouI70aEBtqy7H3f/mstovlMnUvLd5K+9dxH+qUvIV1eb5r9fFyX1AOArL12Fulb2yER8cpgyz+AKhdGbadm69HnHGASMKmVV6JNaKYeGwBvITWl9v0P3J7tu5TXBy20S7O/i7tJX9KKHeEi9Ym1W7P1GA/c2RdcACLiZeBO0g0bW26NDDCNiIg/kTqzX1tEU9Ja2fOr6hz6M1J7Z377AWnM+R7ADysqck2DqAeSNgSuzJ7uHxHLqypnTmULo7ZZs/VA0gmkJphPRVoGaSg1U49lrP4/8I/ZvhOBT1ZT1H41+7v4dfa4R19C1jw9gdcDZjs1W48HgHdlTa4AKN1q/l3A4grK6XkwZTZen6B4EumP7EIKExRJwwNXAbvXOc90OmOiZUP1AK4hDbE+FNglv1VY1lHAX4DZpLXlppI+sL5WyLcI+F4hbRbwJ1IH5gGkEUBDOdGy4Xpk73UAFxTfc2CLbqlHjfOMZWgnWjb7N3V5duz/Jk1QnEua1zOqW+pBCogrSc3c+wH7k4LVSuBvKilru9+cbt1IM6oXAS8BC4A9C/snZv84E+ucYzpDGGCarUf2vOZWcVm3J806XpH9Q50KrF3Is5g0kCKftkn2wfws6Yrxv6kxcbGN73nD9SAF//7e98O6pR41zjFkAWaQf1MbAv8XeDo79lrg3d30N5Wl7UlaIeKZbJtb7zNrsJsXuzQzs0q4D8bMzCrhAGNmZpVwgDEzs0o4wJiZWSUcYMzMrBIOMGZmVgkHGDMzq4QDjJmZVeL/A4S+RQwrpDrYAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# plot estimated values for \"a\"\n",
    "plt.bar(result['values'], result['probabilities'], width=0.5/len(result['probabilities']))\n",
    "plt.xticks([0, 0.25, 0.5, 0.75, 1], size=15)\n",
    "plt.yticks([0, 0.25, 0.5, 0.75, 1], size=15)\n",
    "plt.title('\"a\" Value', size=15)\n",
    "plt.ylabel('Probability', size=15)\n",
    "plt.ylim((0,1))\n",
    "plt.grid()\n",
    "plt.show()\n",
    "\n",
    "# plot estimated values for option price (after re-scaling and reversing the c_approx-transformation)\n",
    "plt.bar(result['mapped_values'], result['probabilities'], width=1/len(result['probabilities']))\n",
    "plt.plot([exact_value, exact_value], [0,1], 'r--', linewidth=2)\n",
    "plt.xticks(size=15)\n",
    "plt.yticks([0, 0.25, 0.5, 0.75, 1], size=15)\n",
    "plt.title('Estimated Option Price', size=15)\n",
    "plt.ylabel('Probability', size=15)\n",
    "plt.ylim((0,1))\n",
    "plt.grid()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}
